1 SUBROUTINE DSYGS2( 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 * DSYGS2 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 * Specifies whether the upper or lower triangular part of the
39 * symmetric matrix A is stored, and how B has been factorized.
40 * = 'U': Upper triangular
41 * = 'L': Lower triangular
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
81 DOUBLE PRECISION AKK, BKK, CT
82 * ..
83 * .. External Subroutines ..
84 EXTERNAL DAXPY, DSCAL, DSYR2, DTRMV, DTRSV, XERBLA
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC MAX
88 * ..
89 * .. External Functions ..
90 LOGICAL LSAME
91 EXTERNAL LSAME
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( 'DSYGS2', -INFO )
112 RETURN
113 END IF
114 *
115 IF( ITYPE.EQ.1 ) THEN
116 IF( UPPER ) THEN
117 *
118 * Compute inv(U**T)*A*inv(U)
119 *
120 DO 10 K = 1, N
121 *
122 * Update the upper triangle of A(k:n,k:n)
123 *
124 AKK = A( K, K )
125 BKK = B( K, K )
126 AKK = AKK / BKK**2
127 A( K, K ) = AKK
128 IF( K.LT.N ) THEN
129 CALL DSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
130 CT = -HALF*AKK
131 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
132 $ LDA )
133 CALL DSYR2( UPLO, N-K, -ONE, A( K, K+1 ), LDA,
134 $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
135 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
136 $ LDA )
137 CALL DTRSV( UPLO, 'Transpose', 'Non-unit', N-K,
138 $ B( K+1, K+1 ), LDB, A( K, K+1 ), LDA )
139 END IF
140 10 CONTINUE
141 ELSE
142 *
143 * Compute inv(L)*A*inv(L**T)
144 *
145 DO 20 K = 1, N
146 *
147 * Update the lower triangle of A(k:n,k:n)
148 *
149 AKK = A( K, K )
150 BKK = B( K, K )
151 AKK = AKK / BKK**2
152 A( K, K ) = AKK
153 IF( K.LT.N ) THEN
154 CALL DSCAL( N-K, ONE / BKK, A( K+1, K ), 1 )
155 CT = -HALF*AKK
156 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
157 CALL DSYR2( UPLO, N-K, -ONE, A( K+1, K ), 1,
158 $ B( K+1, K ), 1, A( K+1, K+1 ), LDA )
159 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
160 CALL DTRSV( UPLO, 'No transpose', 'Non-unit', N-K,
161 $ B( K+1, K+1 ), LDB, A( K+1, K ), 1 )
162 END IF
163 20 CONTINUE
164 END IF
165 ELSE
166 IF( UPPER ) THEN
167 *
168 * Compute U*A*U**T
169 *
170 DO 30 K = 1, N
171 *
172 * Update the upper triangle of A(1:k,1:k)
173 *
174 AKK = A( K, K )
175 BKK = B( K, K )
176 CALL DTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B,
177 $ LDB, A( 1, K ), 1 )
178 CT = HALF*AKK
179 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
180 CALL DSYR2( UPLO, K-1, ONE, A( 1, K ), 1, B( 1, K ), 1,
181 $ A, LDA )
182 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
183 CALL DSCAL( K-1, BKK, A( 1, K ), 1 )
184 A( K, K ) = AKK*BKK**2
185 30 CONTINUE
186 ELSE
187 *
188 * Compute L**T *A*L
189 *
190 DO 40 K = 1, N
191 *
192 * Update the lower triangle of A(1:k,1:k)
193 *
194 AKK = A( K, K )
195 BKK = B( K, K )
196 CALL DTRMV( UPLO, 'Transpose', 'Non-unit', K-1, B, LDB,
197 $ A( K, 1 ), LDA )
198 CT = HALF*AKK
199 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
200 CALL DSYR2( UPLO, K-1, ONE, A( K, 1 ), LDA, B( K, 1 ),
201 $ LDB, A, LDA )
202 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
203 CALL DSCAL( K-1, BKK, A( K, 1 ), LDA )
204 A( K, K ) = AKK*BKK**2
205 40 CONTINUE
206 END IF
207 END IF
208 RETURN
209 *
210 * End of DSYGS2
211 *
212 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 * DSYGS2 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 * Specifies whether the upper or lower triangular part of the
39 * symmetric matrix A is stored, and how B has been factorized.
40 * = 'U': Upper triangular
41 * = 'L': Lower triangular
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
81 DOUBLE PRECISION AKK, BKK, CT
82 * ..
83 * .. External Subroutines ..
84 EXTERNAL DAXPY, DSCAL, DSYR2, DTRMV, DTRSV, XERBLA
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC MAX
88 * ..
89 * .. External Functions ..
90 LOGICAL LSAME
91 EXTERNAL LSAME
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( 'DSYGS2', -INFO )
112 RETURN
113 END IF
114 *
115 IF( ITYPE.EQ.1 ) THEN
116 IF( UPPER ) THEN
117 *
118 * Compute inv(U**T)*A*inv(U)
119 *
120 DO 10 K = 1, N
121 *
122 * Update the upper triangle of A(k:n,k:n)
123 *
124 AKK = A( K, K )
125 BKK = B( K, K )
126 AKK = AKK / BKK**2
127 A( K, K ) = AKK
128 IF( K.LT.N ) THEN
129 CALL DSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
130 CT = -HALF*AKK
131 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
132 $ LDA )
133 CALL DSYR2( UPLO, N-K, -ONE, A( K, K+1 ), LDA,
134 $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
135 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
136 $ LDA )
137 CALL DTRSV( UPLO, 'Transpose', 'Non-unit', N-K,
138 $ B( K+1, K+1 ), LDB, A( K, K+1 ), LDA )
139 END IF
140 10 CONTINUE
141 ELSE
142 *
143 * Compute inv(L)*A*inv(L**T)
144 *
145 DO 20 K = 1, N
146 *
147 * Update the lower triangle of A(k:n,k:n)
148 *
149 AKK = A( K, K )
150 BKK = B( K, K )
151 AKK = AKK / BKK**2
152 A( K, K ) = AKK
153 IF( K.LT.N ) THEN
154 CALL DSCAL( N-K, ONE / BKK, A( K+1, K ), 1 )
155 CT = -HALF*AKK
156 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
157 CALL DSYR2( UPLO, N-K, -ONE, A( K+1, K ), 1,
158 $ B( K+1, K ), 1, A( K+1, K+1 ), LDA )
159 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
160 CALL DTRSV( UPLO, 'No transpose', 'Non-unit', N-K,
161 $ B( K+1, K+1 ), LDB, A( K+1, K ), 1 )
162 END IF
163 20 CONTINUE
164 END IF
165 ELSE
166 IF( UPPER ) THEN
167 *
168 * Compute U*A*U**T
169 *
170 DO 30 K = 1, N
171 *
172 * Update the upper triangle of A(1:k,1:k)
173 *
174 AKK = A( K, K )
175 BKK = B( K, K )
176 CALL DTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B,
177 $ LDB, A( 1, K ), 1 )
178 CT = HALF*AKK
179 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
180 CALL DSYR2( UPLO, K-1, ONE, A( 1, K ), 1, B( 1, K ), 1,
181 $ A, LDA )
182 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
183 CALL DSCAL( K-1, BKK, A( 1, K ), 1 )
184 A( K, K ) = AKK*BKK**2
185 30 CONTINUE
186 ELSE
187 *
188 * Compute L**T *A*L
189 *
190 DO 40 K = 1, N
191 *
192 * Update the lower triangle of A(1:k,1:k)
193 *
194 AKK = A( K, K )
195 BKK = B( K, K )
196 CALL DTRMV( UPLO, 'Transpose', 'Non-unit', K-1, B, LDB,
197 $ A( K, 1 ), LDA )
198 CT = HALF*AKK
199 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
200 CALL DSYR2( UPLO, K-1, ONE, A( K, 1 ), LDA, B( K, 1 ),
201 $ LDB, A, LDA )
202 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
203 CALL DSCAL( K-1, BKK, A( K, 1 ), LDA )
204 A( K, K ) = AKK*BKK**2
205 40 CONTINUE
206 END IF
207 END IF
208 RETURN
209 *
210 * End of DSYGS2
211 *
212 END