1 SUBROUTINE DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, 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, LDA, LDB, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), B( LDB, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSYGST reduces a real symmetric-definite generalized eigenproblem
20 * to standard form.
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 DPOTRF.
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 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
47 * On entry, the symmetric matrix A. If UPLO = 'U', the leading
48 * N-by-N upper triangular part of A contains the upper
49 * triangular part of the matrix A, and the strictly lower
50 * triangular part of A is not referenced. If UPLO = 'L', the
51 * leading N-by-N lower triangular part of A contains the lower
52 * triangular part of the matrix A, and the strictly upper
53 * triangular part of A is not referenced.
54 *
55 * On exit, if INFO = 0, the transformed matrix, stored in the
56 * same format as A.
57 *
58 * LDA (input) INTEGER
59 * The leading dimension of the array A. LDA >= max(1,N).
60 *
61 * B (input) DOUBLE PRECISION array, dimension (LDB,N)
62 * The triangular factor from the Cholesky factorization of B,
63 * as returned by DPOTRF.
64 *
65 * LDB (input) INTEGER
66 * The leading dimension of the array B. LDB >= max(1,N).
67 *
68 * INFO (output) INTEGER
69 * = 0: successful exit
70 * < 0: if INFO = -i, the i-th argument had an illegal value
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 DOUBLE PRECISION ONE, HALF
76 PARAMETER ( ONE = 1.0D0, HALF = 0.5D0 )
77 * ..
78 * .. Local Scalars ..
79 LOGICAL UPPER
80 INTEGER K, KB, NB
81 * ..
82 * .. External Subroutines ..
83 EXTERNAL DSYGS2, DSYMM, DSYR2K, DTRMM, DTRSM, XERBLA
84 * ..
85 * .. Intrinsic Functions ..
86 INTRINSIC MAX, MIN
87 * ..
88 * .. External Functions ..
89 LOGICAL LSAME
90 INTEGER ILAENV
91 EXTERNAL LSAME, ILAENV
92 * ..
93 * .. Executable Statements ..
94 *
95 * Test the input parameters.
96 *
97 INFO = 0
98 UPPER = LSAME( UPLO, 'U' )
99 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
100 INFO = -1
101 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
102 INFO = -2
103 ELSE IF( N.LT.0 ) THEN
104 INFO = -3
105 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
106 INFO = -5
107 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
108 INFO = -7
109 END IF
110 IF( INFO.NE.0 ) THEN
111 CALL XERBLA( 'DSYGST', -INFO )
112 RETURN
113 END IF
114 *
115 * Quick return if possible
116 *
117 IF( N.EQ.0 )
118 $ RETURN
119 *
120 * Determine the block size for this environment.
121 *
122 NB = ILAENV( 1, 'DSYGST', UPLO, N, -1, -1, -1 )
123 *
124 IF( NB.LE.1 .OR. NB.GE.N ) THEN
125 *
126 * Use unblocked code
127 *
128 CALL DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
129 ELSE
130 *
131 * Use blocked code
132 *
133 IF( ITYPE.EQ.1 ) THEN
134 IF( UPPER ) THEN
135 *
136 * Compute inv(U**T)*A*inv(U)
137 *
138 DO 10 K = 1, N, NB
139 KB = MIN( N-K+1, NB )
140 *
141 * Update the upper triangle of A(k:n,k:n)
142 *
143 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
144 $ B( K, K ), LDB, INFO )
145 IF( K+KB.LE.N ) THEN
146 CALL DTRSM( 'Left', UPLO, 'Transpose', 'Non-unit',
147 $ KB, N-K-KB+1, ONE, B( K, K ), LDB,
148 $ A( K, K+KB ), LDA )
149 CALL DSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
150 $ A( K, K ), LDA, B( K, K+KB ), LDB, ONE,
151 $ A( K, K+KB ), LDA )
152 CALL DSYR2K( UPLO, 'Transpose', N-K-KB+1, KB, -ONE,
153 $ A( K, K+KB ), LDA, B( K, K+KB ), LDB,
154 $ ONE, A( K+KB, K+KB ), LDA )
155 CALL DSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
156 $ A( K, K ), LDA, B( K, K+KB ), LDB, ONE,
157 $ A( K, K+KB ), LDA )
158 CALL DTRSM( 'Right', UPLO, 'No transpose',
159 $ 'Non-unit', KB, N-K-KB+1, ONE,
160 $ B( K+KB, K+KB ), LDB, A( K, K+KB ),
161 $ LDA )
162 END IF
163 10 CONTINUE
164 ELSE
165 *
166 * Compute inv(L)*A*inv(L**T)
167 *
168 DO 20 K = 1, N, NB
169 KB = MIN( N-K+1, NB )
170 *
171 * Update the lower triangle of A(k:n,k:n)
172 *
173 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
174 $ B( K, K ), LDB, INFO )
175 IF( K+KB.LE.N ) THEN
176 CALL DTRSM( 'Right', UPLO, 'Transpose', 'Non-unit',
177 $ N-K-KB+1, KB, ONE, B( K, K ), LDB,
178 $ A( K+KB, K ), LDA )
179 CALL DSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
180 $ A( K, K ), LDA, B( K+KB, K ), LDB, ONE,
181 $ A( K+KB, K ), LDA )
182 CALL DSYR2K( UPLO, 'No transpose', N-K-KB+1, KB,
183 $ -ONE, A( K+KB, K ), LDA, B( K+KB, K ),
184 $ LDB, ONE, A( K+KB, K+KB ), LDA )
185 CALL DSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
186 $ A( K, K ), LDA, B( K+KB, K ), LDB, ONE,
187 $ A( K+KB, K ), LDA )
188 CALL DTRSM( 'Left', UPLO, 'No transpose',
189 $ 'Non-unit', N-K-KB+1, KB, ONE,
190 $ B( K+KB, K+KB ), LDB, A( K+KB, K ),
191 $ LDA )
192 END IF
193 20 CONTINUE
194 END IF
195 ELSE
196 IF( UPPER ) THEN
197 *
198 * Compute U*A*U**T
199 *
200 DO 30 K = 1, N, NB
201 KB = MIN( N-K+1, NB )
202 *
203 * Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
204 *
205 CALL DTRMM( 'Left', UPLO, 'No transpose', 'Non-unit',
206 $ K-1, KB, ONE, B, LDB, A( 1, K ), LDA )
207 CALL DSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
208 $ LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA )
209 CALL DSYR2K( UPLO, 'No transpose', K-1, KB, ONE,
210 $ A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
211 $ LDA )
212 CALL DSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
213 $ LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA )
214 CALL DTRMM( 'Right', UPLO, 'Transpose', 'Non-unit',
215 $ K-1, KB, ONE, B( K, K ), LDB, A( 1, K ),
216 $ LDA )
217 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
218 $ B( K, K ), LDB, INFO )
219 30 CONTINUE
220 ELSE
221 *
222 * Compute L**T*A*L
223 *
224 DO 40 K = 1, N, NB
225 KB = MIN( N-K+1, NB )
226 *
227 * Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
228 *
229 CALL DTRMM( 'Right', UPLO, 'No transpose', 'Non-unit',
230 $ KB, K-1, ONE, B, LDB, A( K, 1 ), LDA )
231 CALL DSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
232 $ LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA )
233 CALL DSYR2K( UPLO, 'Transpose', K-1, KB, ONE,
234 $ A( K, 1 ), LDA, B( K, 1 ), LDB, ONE, A,
235 $ LDA )
236 CALL DSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
237 $ LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA )
238 CALL DTRMM( 'Left', UPLO, 'Transpose', 'Non-unit', KB,
239 $ K-1, ONE, B( K, K ), LDB, A( K, 1 ), LDA )
240 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
241 $ B( K, K ), LDB, INFO )
242 40 CONTINUE
243 END IF
244 END IF
245 END IF
246 RETURN
247 *
248 * End of DSYGST
249 *
250 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, LDA, LDB, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), B( LDB, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSYGST reduces a real symmetric-definite generalized eigenproblem
20 * to standard form.
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 DPOTRF.
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 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
47 * On entry, the symmetric matrix A. If UPLO = 'U', the leading
48 * N-by-N upper triangular part of A contains the upper
49 * triangular part of the matrix A, and the strictly lower
50 * triangular part of A is not referenced. If UPLO = 'L', the
51 * leading N-by-N lower triangular part of A contains the lower
52 * triangular part of the matrix A, and the strictly upper
53 * triangular part of A is not referenced.
54 *
55 * On exit, if INFO = 0, the transformed matrix, stored in the
56 * same format as A.
57 *
58 * LDA (input) INTEGER
59 * The leading dimension of the array A. LDA >= max(1,N).
60 *
61 * B (input) DOUBLE PRECISION array, dimension (LDB,N)
62 * The triangular factor from the Cholesky factorization of B,
63 * as returned by DPOTRF.
64 *
65 * LDB (input) INTEGER
66 * The leading dimension of the array B. LDB >= max(1,N).
67 *
68 * INFO (output) INTEGER
69 * = 0: successful exit
70 * < 0: if INFO = -i, the i-th argument had an illegal value
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 DOUBLE PRECISION ONE, HALF
76 PARAMETER ( ONE = 1.0D0, HALF = 0.5D0 )
77 * ..
78 * .. Local Scalars ..
79 LOGICAL UPPER
80 INTEGER K, KB, NB
81 * ..
82 * .. External Subroutines ..
83 EXTERNAL DSYGS2, DSYMM, DSYR2K, DTRMM, DTRSM, XERBLA
84 * ..
85 * .. Intrinsic Functions ..
86 INTRINSIC MAX, MIN
87 * ..
88 * .. External Functions ..
89 LOGICAL LSAME
90 INTEGER ILAENV
91 EXTERNAL LSAME, ILAENV
92 * ..
93 * .. Executable Statements ..
94 *
95 * Test the input parameters.
96 *
97 INFO = 0
98 UPPER = LSAME( UPLO, 'U' )
99 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
100 INFO = -1
101 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
102 INFO = -2
103 ELSE IF( N.LT.0 ) THEN
104 INFO = -3
105 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
106 INFO = -5
107 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
108 INFO = -7
109 END IF
110 IF( INFO.NE.0 ) THEN
111 CALL XERBLA( 'DSYGST', -INFO )
112 RETURN
113 END IF
114 *
115 * Quick return if possible
116 *
117 IF( N.EQ.0 )
118 $ RETURN
119 *
120 * Determine the block size for this environment.
121 *
122 NB = ILAENV( 1, 'DSYGST', UPLO, N, -1, -1, -1 )
123 *
124 IF( NB.LE.1 .OR. NB.GE.N ) THEN
125 *
126 * Use unblocked code
127 *
128 CALL DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
129 ELSE
130 *
131 * Use blocked code
132 *
133 IF( ITYPE.EQ.1 ) THEN
134 IF( UPPER ) THEN
135 *
136 * Compute inv(U**T)*A*inv(U)
137 *
138 DO 10 K = 1, N, NB
139 KB = MIN( N-K+1, NB )
140 *
141 * Update the upper triangle of A(k:n,k:n)
142 *
143 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
144 $ B( K, K ), LDB, INFO )
145 IF( K+KB.LE.N ) THEN
146 CALL DTRSM( 'Left', UPLO, 'Transpose', 'Non-unit',
147 $ KB, N-K-KB+1, ONE, B( K, K ), LDB,
148 $ A( K, K+KB ), LDA )
149 CALL DSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
150 $ A( K, K ), LDA, B( K, K+KB ), LDB, ONE,
151 $ A( K, K+KB ), LDA )
152 CALL DSYR2K( UPLO, 'Transpose', N-K-KB+1, KB, -ONE,
153 $ A( K, K+KB ), LDA, B( K, K+KB ), LDB,
154 $ ONE, A( K+KB, K+KB ), LDA )
155 CALL DSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
156 $ A( K, K ), LDA, B( K, K+KB ), LDB, ONE,
157 $ A( K, K+KB ), LDA )
158 CALL DTRSM( 'Right', UPLO, 'No transpose',
159 $ 'Non-unit', KB, N-K-KB+1, ONE,
160 $ B( K+KB, K+KB ), LDB, A( K, K+KB ),
161 $ LDA )
162 END IF
163 10 CONTINUE
164 ELSE
165 *
166 * Compute inv(L)*A*inv(L**T)
167 *
168 DO 20 K = 1, N, NB
169 KB = MIN( N-K+1, NB )
170 *
171 * Update the lower triangle of A(k:n,k:n)
172 *
173 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
174 $ B( K, K ), LDB, INFO )
175 IF( K+KB.LE.N ) THEN
176 CALL DTRSM( 'Right', UPLO, 'Transpose', 'Non-unit',
177 $ N-K-KB+1, KB, ONE, B( K, K ), LDB,
178 $ A( K+KB, K ), LDA )
179 CALL DSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
180 $ A( K, K ), LDA, B( K+KB, K ), LDB, ONE,
181 $ A( K+KB, K ), LDA )
182 CALL DSYR2K( UPLO, 'No transpose', N-K-KB+1, KB,
183 $ -ONE, A( K+KB, K ), LDA, B( K+KB, K ),
184 $ LDB, ONE, A( K+KB, K+KB ), LDA )
185 CALL DSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
186 $ A( K, K ), LDA, B( K+KB, K ), LDB, ONE,
187 $ A( K+KB, K ), LDA )
188 CALL DTRSM( 'Left', UPLO, 'No transpose',
189 $ 'Non-unit', N-K-KB+1, KB, ONE,
190 $ B( K+KB, K+KB ), LDB, A( K+KB, K ),
191 $ LDA )
192 END IF
193 20 CONTINUE
194 END IF
195 ELSE
196 IF( UPPER ) THEN
197 *
198 * Compute U*A*U**T
199 *
200 DO 30 K = 1, N, NB
201 KB = MIN( N-K+1, NB )
202 *
203 * Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
204 *
205 CALL DTRMM( 'Left', UPLO, 'No transpose', 'Non-unit',
206 $ K-1, KB, ONE, B, LDB, A( 1, K ), LDA )
207 CALL DSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
208 $ LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA )
209 CALL DSYR2K( UPLO, 'No transpose', K-1, KB, ONE,
210 $ A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
211 $ LDA )
212 CALL DSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
213 $ LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA )
214 CALL DTRMM( 'Right', UPLO, 'Transpose', 'Non-unit',
215 $ K-1, KB, ONE, B( K, K ), LDB, A( 1, K ),
216 $ LDA )
217 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
218 $ B( K, K ), LDB, INFO )
219 30 CONTINUE
220 ELSE
221 *
222 * Compute L**T*A*L
223 *
224 DO 40 K = 1, N, NB
225 KB = MIN( N-K+1, NB )
226 *
227 * Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
228 *
229 CALL DTRMM( 'Right', UPLO, 'No transpose', 'Non-unit',
230 $ KB, K-1, ONE, B, LDB, A( K, 1 ), LDA )
231 CALL DSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
232 $ LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA )
233 CALL DSYR2K( UPLO, 'Transpose', K-1, KB, ONE,
234 $ A( K, 1 ), LDA, B( K, 1 ), LDB, ONE, A,
235 $ LDA )
236 CALL DSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
237 $ LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA )
238 CALL DTRMM( 'Left', UPLO, 'Transpose', 'Non-unit', KB,
239 $ K-1, ONE, B( K, K ), LDB, A( K, 1 ), LDA )
240 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
241 $ B( K, K ), LDB, INFO )
242 40 CONTINUE
243 END IF
244 END IF
245 END IF
246 RETURN
247 *
248 * End of DSYGST
249 *
250 END