1 SUBROUTINE DSPGST( ITYPE, UPLO, N, AP, BP, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, ITYPE, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * ), BP( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSPGST reduces a real symmetric-definite generalized eigenproblem
20 * to standard form, using packed storage.
21 *
22 * If ITYPE = 1, the problem is A*x = lambda*B*x,
23 * and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
24 *
25 * If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
26 * B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
27 *
28 * B must have been previously factorized as U**T*U or L*L**T by DPPTRF.
29 *
30 * Arguments
31 * =========
32 *
33 * ITYPE (input) INTEGER
34 * = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
35 * = 2 or 3: compute U*A*U**T or L**T*A*L.
36 *
37 * UPLO (input) CHARACTER*1
38 * = 'U': Upper triangle of A is stored and B is factored as
39 * U**T*U;
40 * = 'L': Lower triangle of A is stored and B is factored as
41 * L*L**T.
42 *
43 * N (input) INTEGER
44 * The order of the matrices A and B. N >= 0.
45 *
46 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
47 * On entry, the upper or lower triangle of the symmetric matrix
48 * A, packed columnwise in a linear array. The j-th column of A
49 * is stored in the array AP as follows:
50 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
51 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
52 *
53 * On exit, if INFO = 0, the transformed matrix, stored in the
54 * same format as A.
55 *
56 * BP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
57 * The triangular factor from the Cholesky factorization of B,
58 * stored in the same format as A, as returned by DPPTRF.
59 *
60 * INFO (output) INTEGER
61 * = 0: successful exit
62 * < 0: if INFO = -i, the i-th argument had an illegal value
63 *
64 * =====================================================================
65 *
66 * .. Parameters ..
67 DOUBLE PRECISION ONE, HALF
68 PARAMETER ( ONE = 1.0D0, HALF = 0.5D0 )
69 * ..
70 * .. Local Scalars ..
71 LOGICAL UPPER
72 INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK
73 DOUBLE PRECISION AJJ, AKK, BJJ, BKK, CT
74 * ..
75 * .. External Subroutines ..
76 EXTERNAL DAXPY, DSCAL, DSPMV, DSPR2, DTPMV, DTPSV,
77 $ XERBLA
78 * ..
79 * .. External Functions ..
80 LOGICAL LSAME
81 DOUBLE PRECISION DDOT
82 EXTERNAL LSAME, DDOT
83 * ..
84 * .. Executable Statements ..
85 *
86 * Test the input parameters.
87 *
88 INFO = 0
89 UPPER = LSAME( UPLO, 'U' )
90 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
91 INFO = -1
92 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
93 INFO = -2
94 ELSE IF( N.LT.0 ) THEN
95 INFO = -3
96 END IF
97 IF( INFO.NE.0 ) THEN
98 CALL XERBLA( 'DSPGST', -INFO )
99 RETURN
100 END IF
101 *
102 IF( ITYPE.EQ.1 ) THEN
103 IF( UPPER ) THEN
104 *
105 * Compute inv(U**T)*A*inv(U)
106 *
107 * J1 and JJ are the indices of A(1,j) and A(j,j)
108 *
109 JJ = 0
110 DO 10 J = 1, N
111 J1 = JJ + 1
112 JJ = JJ + J
113 *
114 * Compute the j-th column of the upper triangle of A
115 *
116 BJJ = BP( JJ )
117 CALL DTPSV( UPLO, 'Transpose', 'Nonunit', J, BP,
118 $ AP( J1 ), 1 )
119 CALL DSPMV( UPLO, J-1, -ONE, AP, BP( J1 ), 1, ONE,
120 $ AP( J1 ), 1 )
121 CALL DSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
122 AP( JJ ) = ( AP( JJ )-DDOT( J-1, AP( J1 ), 1, BP( J1 ),
123 $ 1 ) ) / BJJ
124 10 CONTINUE
125 ELSE
126 *
127 * Compute inv(L)*A*inv(L**T)
128 *
129 * KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
130 *
131 KK = 1
132 DO 20 K = 1, N
133 K1K1 = KK + N - K + 1
134 *
135 * Update the lower triangle of A(k:n,k:n)
136 *
137 AKK = AP( KK )
138 BKK = BP( KK )
139 AKK = AKK / BKK**2
140 AP( KK ) = AKK
141 IF( K.LT.N ) THEN
142 CALL DSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
143 CT = -HALF*AKK
144 CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
145 CALL DSPR2( UPLO, N-K, -ONE, AP( KK+1 ), 1,
146 $ BP( KK+1 ), 1, AP( K1K1 ) )
147 CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
148 CALL DTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
149 $ BP( K1K1 ), AP( KK+1 ), 1 )
150 END IF
151 KK = K1K1
152 20 CONTINUE
153 END IF
154 ELSE
155 IF( UPPER ) THEN
156 *
157 * Compute U*A*U**T
158 *
159 * K1 and KK are the indices of A(1,k) and A(k,k)
160 *
161 KK = 0
162 DO 30 K = 1, N
163 K1 = KK + 1
164 KK = KK + K
165 *
166 * Update the upper triangle of A(1:k,1:k)
167 *
168 AKK = AP( KK )
169 BKK = BP( KK )
170 CALL DTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
171 $ AP( K1 ), 1 )
172 CT = HALF*AKK
173 CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
174 CALL DSPR2( UPLO, K-1, ONE, AP( K1 ), 1, BP( K1 ), 1,
175 $ AP )
176 CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
177 CALL DSCAL( K-1, BKK, AP( K1 ), 1 )
178 AP( KK ) = AKK*BKK**2
179 30 CONTINUE
180 ELSE
181 *
182 * Compute L**T *A*L
183 *
184 * JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
185 *
186 JJ = 1
187 DO 40 J = 1, N
188 J1J1 = JJ + N - J + 1
189 *
190 * Compute the j-th column of the lower triangle of A
191 *
192 AJJ = AP( JJ )
193 BJJ = BP( JJ )
194 AP( JJ ) = AJJ*BJJ + DDOT( N-J, AP( JJ+1 ), 1,
195 $ BP( JJ+1 ), 1 )
196 CALL DSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
197 CALL DSPMV( UPLO, N-J, ONE, AP( J1J1 ), BP( JJ+1 ), 1,
198 $ ONE, AP( JJ+1 ), 1 )
199 CALL DTPMV( UPLO, 'Transpose', 'Non-unit', N-J+1,
200 $ BP( JJ ), AP( JJ ), 1 )
201 JJ = J1J1
202 40 CONTINUE
203 END IF
204 END IF
205 RETURN
206 *
207 * End of DSPGST
208 *
209 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, ITYPE, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * ), BP( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSPGST reduces a real symmetric-definite generalized eigenproblem
20 * to standard form, using packed storage.
21 *
22 * If ITYPE = 1, the problem is A*x = lambda*B*x,
23 * and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
24 *
25 * If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
26 * B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
27 *
28 * B must have been previously factorized as U**T*U or L*L**T by DPPTRF.
29 *
30 * Arguments
31 * =========
32 *
33 * ITYPE (input) INTEGER
34 * = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
35 * = 2 or 3: compute U*A*U**T or L**T*A*L.
36 *
37 * UPLO (input) CHARACTER*1
38 * = 'U': Upper triangle of A is stored and B is factored as
39 * U**T*U;
40 * = 'L': Lower triangle of A is stored and B is factored as
41 * L*L**T.
42 *
43 * N (input) INTEGER
44 * The order of the matrices A and B. N >= 0.
45 *
46 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
47 * On entry, the upper or lower triangle of the symmetric matrix
48 * A, packed columnwise in a linear array. The j-th column of A
49 * is stored in the array AP as follows:
50 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
51 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
52 *
53 * On exit, if INFO = 0, the transformed matrix, stored in the
54 * same format as A.
55 *
56 * BP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
57 * The triangular factor from the Cholesky factorization of B,
58 * stored in the same format as A, as returned by DPPTRF.
59 *
60 * INFO (output) INTEGER
61 * = 0: successful exit
62 * < 0: if INFO = -i, the i-th argument had an illegal value
63 *
64 * =====================================================================
65 *
66 * .. Parameters ..
67 DOUBLE PRECISION ONE, HALF
68 PARAMETER ( ONE = 1.0D0, HALF = 0.5D0 )
69 * ..
70 * .. Local Scalars ..
71 LOGICAL UPPER
72 INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK
73 DOUBLE PRECISION AJJ, AKK, BJJ, BKK, CT
74 * ..
75 * .. External Subroutines ..
76 EXTERNAL DAXPY, DSCAL, DSPMV, DSPR2, DTPMV, DTPSV,
77 $ XERBLA
78 * ..
79 * .. External Functions ..
80 LOGICAL LSAME
81 DOUBLE PRECISION DDOT
82 EXTERNAL LSAME, DDOT
83 * ..
84 * .. Executable Statements ..
85 *
86 * Test the input parameters.
87 *
88 INFO = 0
89 UPPER = LSAME( UPLO, 'U' )
90 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
91 INFO = -1
92 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
93 INFO = -2
94 ELSE IF( N.LT.0 ) THEN
95 INFO = -3
96 END IF
97 IF( INFO.NE.0 ) THEN
98 CALL XERBLA( 'DSPGST', -INFO )
99 RETURN
100 END IF
101 *
102 IF( ITYPE.EQ.1 ) THEN
103 IF( UPPER ) THEN
104 *
105 * Compute inv(U**T)*A*inv(U)
106 *
107 * J1 and JJ are the indices of A(1,j) and A(j,j)
108 *
109 JJ = 0
110 DO 10 J = 1, N
111 J1 = JJ + 1
112 JJ = JJ + J
113 *
114 * Compute the j-th column of the upper triangle of A
115 *
116 BJJ = BP( JJ )
117 CALL DTPSV( UPLO, 'Transpose', 'Nonunit', J, BP,
118 $ AP( J1 ), 1 )
119 CALL DSPMV( UPLO, J-1, -ONE, AP, BP( J1 ), 1, ONE,
120 $ AP( J1 ), 1 )
121 CALL DSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
122 AP( JJ ) = ( AP( JJ )-DDOT( J-1, AP( J1 ), 1, BP( J1 ),
123 $ 1 ) ) / BJJ
124 10 CONTINUE
125 ELSE
126 *
127 * Compute inv(L)*A*inv(L**T)
128 *
129 * KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
130 *
131 KK = 1
132 DO 20 K = 1, N
133 K1K1 = KK + N - K + 1
134 *
135 * Update the lower triangle of A(k:n,k:n)
136 *
137 AKK = AP( KK )
138 BKK = BP( KK )
139 AKK = AKK / BKK**2
140 AP( KK ) = AKK
141 IF( K.LT.N ) THEN
142 CALL DSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
143 CT = -HALF*AKK
144 CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
145 CALL DSPR2( UPLO, N-K, -ONE, AP( KK+1 ), 1,
146 $ BP( KK+1 ), 1, AP( K1K1 ) )
147 CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
148 CALL DTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
149 $ BP( K1K1 ), AP( KK+1 ), 1 )
150 END IF
151 KK = K1K1
152 20 CONTINUE
153 END IF
154 ELSE
155 IF( UPPER ) THEN
156 *
157 * Compute U*A*U**T
158 *
159 * K1 and KK are the indices of A(1,k) and A(k,k)
160 *
161 KK = 0
162 DO 30 K = 1, N
163 K1 = KK + 1
164 KK = KK + K
165 *
166 * Update the upper triangle of A(1:k,1:k)
167 *
168 AKK = AP( KK )
169 BKK = BP( KK )
170 CALL DTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
171 $ AP( K1 ), 1 )
172 CT = HALF*AKK
173 CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
174 CALL DSPR2( UPLO, K-1, ONE, AP( K1 ), 1, BP( K1 ), 1,
175 $ AP )
176 CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
177 CALL DSCAL( K-1, BKK, AP( K1 ), 1 )
178 AP( KK ) = AKK*BKK**2
179 30 CONTINUE
180 ELSE
181 *
182 * Compute L**T *A*L
183 *
184 * JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
185 *
186 JJ = 1
187 DO 40 J = 1, N
188 J1J1 = JJ + N - J + 1
189 *
190 * Compute the j-th column of the lower triangle of A
191 *
192 AJJ = AP( JJ )
193 BJJ = BP( JJ )
194 AP( JJ ) = AJJ*BJJ + DDOT( N-J, AP( JJ+1 ), 1,
195 $ BP( JJ+1 ), 1 )
196 CALL DSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
197 CALL DSPMV( UPLO, N-J, ONE, AP( J1J1 ), BP( JJ+1 ), 1,
198 $ ONE, AP( JJ+1 ), 1 )
199 CALL DTPMV( UPLO, 'Transpose', 'Non-unit', N-J+1,
200 $ BP( JJ ), AP( JJ ), 1 )
201 JJ = J1J1
202 40 CONTINUE
203 END IF
204 END IF
205 RETURN
206 *
207 * End of DSPGST
208 *
209 END