1 SUBROUTINE ZHEGS2( 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 COMPLEX*16 A( LDA, * ), B( LDB, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZHEGS2 reduces a complex Hermitian-definite generalized
20 * eigenproblem to standard form.
21 *
22 * If ITYPE = 1, the problem is A*x = lambda*B*x,
23 * and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
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**H or L**H *A*L.
27 *
28 * B must have been previously factorized as U**H *U or L*L**H by ZPOTRF.
29 *
30 * Arguments
31 * =========
32 *
33 * ITYPE (input) INTEGER
34 * = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
35 * = 2 or 3: compute U*A*U**H or L**H *A*L.
36 *
37 * UPLO (input) CHARACTER*1
38 * Specifies whether the upper or lower triangular part of the
39 * Hermitian 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) COMPLEX*16 array, dimension (LDA,N)
47 * On entry, the Hermitian 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) COMPLEX*16 array, dimension (LDB,N)
62 * The triangular factor from the Cholesky factorization of B,
63 * as returned by ZPOTRF.
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.0D+0, HALF = 0.5D+0 )
77 COMPLEX*16 CONE
78 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
79 * ..
80 * .. Local Scalars ..
81 LOGICAL UPPER
82 INTEGER K
83 DOUBLE PRECISION AKK, BKK
84 COMPLEX*16 CT
85 * ..
86 * .. External Subroutines ..
87 EXTERNAL XERBLA, ZAXPY, ZDSCAL, ZHER2, ZLACGV, ZTRMV,
88 $ ZTRSV
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC MAX
92 * ..
93 * .. External Functions ..
94 LOGICAL LSAME
95 EXTERNAL LSAME
96 * ..
97 * .. Executable Statements ..
98 *
99 * Test the input parameters.
100 *
101 INFO = 0
102 UPPER = LSAME( UPLO, 'U' )
103 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
104 INFO = -1
105 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
106 INFO = -2
107 ELSE IF( N.LT.0 ) THEN
108 INFO = -3
109 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
110 INFO = -5
111 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
112 INFO = -7
113 END IF
114 IF( INFO.NE.0 ) THEN
115 CALL XERBLA( 'ZHEGS2', -INFO )
116 RETURN
117 END IF
118 *
119 IF( ITYPE.EQ.1 ) THEN
120 IF( UPPER ) THEN
121 *
122 * Compute inv(U**H)*A*inv(U)
123 *
124 DO 10 K = 1, N
125 *
126 * Update the upper triangle of A(k:n,k:n)
127 *
128 AKK = A( K, K )
129 BKK = B( K, K )
130 AKK = AKK / BKK**2
131 A( K, K ) = AKK
132 IF( K.LT.N ) THEN
133 CALL ZDSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
134 CT = -HALF*AKK
135 CALL ZLACGV( N-K, A( K, K+1 ), LDA )
136 CALL ZLACGV( N-K, B( K, K+1 ), LDB )
137 CALL ZAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
138 $ LDA )
139 CALL ZHER2( UPLO, N-K, -CONE, A( K, K+1 ), LDA,
140 $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
141 CALL ZAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
142 $ LDA )
143 CALL ZLACGV( N-K, B( K, K+1 ), LDB )
144 CALL ZTRSV( UPLO, 'Conjugate transpose', 'Non-unit',
145 $ N-K, B( K+1, K+1 ), LDB, A( K, K+1 ),
146 $ LDA )
147 CALL ZLACGV( N-K, A( K, K+1 ), LDA )
148 END IF
149 10 CONTINUE
150 ELSE
151 *
152 * Compute inv(L)*A*inv(L**H)
153 *
154 DO 20 K = 1, N
155 *
156 * Update the lower triangle of A(k:n,k:n)
157 *
158 AKK = A( K, K )
159 BKK = B( K, K )
160 AKK = AKK / BKK**2
161 A( K, K ) = AKK
162 IF( K.LT.N ) THEN
163 CALL ZDSCAL( N-K, ONE / BKK, A( K+1, K ), 1 )
164 CT = -HALF*AKK
165 CALL ZAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
166 CALL ZHER2( UPLO, N-K, -CONE, A( K+1, K ), 1,
167 $ B( K+1, K ), 1, A( K+1, K+1 ), LDA )
168 CALL ZAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
169 CALL ZTRSV( UPLO, 'No transpose', 'Non-unit', N-K,
170 $ B( K+1, K+1 ), LDB, A( K+1, K ), 1 )
171 END IF
172 20 CONTINUE
173 END IF
174 ELSE
175 IF( UPPER ) THEN
176 *
177 * Compute U*A*U**H
178 *
179 DO 30 K = 1, N
180 *
181 * Update the upper triangle of A(1:k,1:k)
182 *
183 AKK = A( K, K )
184 BKK = B( K, K )
185 CALL ZTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B,
186 $ LDB, A( 1, K ), 1 )
187 CT = HALF*AKK
188 CALL ZAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
189 CALL ZHER2( UPLO, K-1, CONE, A( 1, K ), 1, B( 1, K ), 1,
190 $ A, LDA )
191 CALL ZAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
192 CALL ZDSCAL( K-1, BKK, A( 1, K ), 1 )
193 A( K, K ) = AKK*BKK**2
194 30 CONTINUE
195 ELSE
196 *
197 * Compute L**H *A*L
198 *
199 DO 40 K = 1, N
200 *
201 * Update the lower triangle of A(1:k,1:k)
202 *
203 AKK = A( K, K )
204 BKK = B( K, K )
205 CALL ZLACGV( K-1, A( K, 1 ), LDA )
206 CALL ZTRMV( UPLO, 'Conjugate transpose', 'Non-unit', K-1,
207 $ B, LDB, A( K, 1 ), LDA )
208 CT = HALF*AKK
209 CALL ZLACGV( K-1, B( K, 1 ), LDB )
210 CALL ZAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
211 CALL ZHER2( UPLO, K-1, CONE, A( K, 1 ), LDA, B( K, 1 ),
212 $ LDB, A, LDA )
213 CALL ZAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
214 CALL ZLACGV( K-1, B( K, 1 ), LDB )
215 CALL ZDSCAL( K-1, BKK, A( K, 1 ), LDA )
216 CALL ZLACGV( K-1, A( K, 1 ), LDA )
217 A( K, K ) = AKK*BKK**2
218 40 CONTINUE
219 END IF
220 END IF
221 RETURN
222 *
223 * End of ZHEGS2
224 *
225 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 COMPLEX*16 A( LDA, * ), B( LDB, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZHEGS2 reduces a complex Hermitian-definite generalized
20 * eigenproblem to standard form.
21 *
22 * If ITYPE = 1, the problem is A*x = lambda*B*x,
23 * and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
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**H or L**H *A*L.
27 *
28 * B must have been previously factorized as U**H *U or L*L**H by ZPOTRF.
29 *
30 * Arguments
31 * =========
32 *
33 * ITYPE (input) INTEGER
34 * = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
35 * = 2 or 3: compute U*A*U**H or L**H *A*L.
36 *
37 * UPLO (input) CHARACTER*1
38 * Specifies whether the upper or lower triangular part of the
39 * Hermitian 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) COMPLEX*16 array, dimension (LDA,N)
47 * On entry, the Hermitian 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) COMPLEX*16 array, dimension (LDB,N)
62 * The triangular factor from the Cholesky factorization of B,
63 * as returned by ZPOTRF.
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.0D+0, HALF = 0.5D+0 )
77 COMPLEX*16 CONE
78 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
79 * ..
80 * .. Local Scalars ..
81 LOGICAL UPPER
82 INTEGER K
83 DOUBLE PRECISION AKK, BKK
84 COMPLEX*16 CT
85 * ..
86 * .. External Subroutines ..
87 EXTERNAL XERBLA, ZAXPY, ZDSCAL, ZHER2, ZLACGV, ZTRMV,
88 $ ZTRSV
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC MAX
92 * ..
93 * .. External Functions ..
94 LOGICAL LSAME
95 EXTERNAL LSAME
96 * ..
97 * .. Executable Statements ..
98 *
99 * Test the input parameters.
100 *
101 INFO = 0
102 UPPER = LSAME( UPLO, 'U' )
103 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
104 INFO = -1
105 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
106 INFO = -2
107 ELSE IF( N.LT.0 ) THEN
108 INFO = -3
109 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
110 INFO = -5
111 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
112 INFO = -7
113 END IF
114 IF( INFO.NE.0 ) THEN
115 CALL XERBLA( 'ZHEGS2', -INFO )
116 RETURN
117 END IF
118 *
119 IF( ITYPE.EQ.1 ) THEN
120 IF( UPPER ) THEN
121 *
122 * Compute inv(U**H)*A*inv(U)
123 *
124 DO 10 K = 1, N
125 *
126 * Update the upper triangle of A(k:n,k:n)
127 *
128 AKK = A( K, K )
129 BKK = B( K, K )
130 AKK = AKK / BKK**2
131 A( K, K ) = AKK
132 IF( K.LT.N ) THEN
133 CALL ZDSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
134 CT = -HALF*AKK
135 CALL ZLACGV( N-K, A( K, K+1 ), LDA )
136 CALL ZLACGV( N-K, B( K, K+1 ), LDB )
137 CALL ZAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
138 $ LDA )
139 CALL ZHER2( UPLO, N-K, -CONE, A( K, K+1 ), LDA,
140 $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
141 CALL ZAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
142 $ LDA )
143 CALL ZLACGV( N-K, B( K, K+1 ), LDB )
144 CALL ZTRSV( UPLO, 'Conjugate transpose', 'Non-unit',
145 $ N-K, B( K+1, K+1 ), LDB, A( K, K+1 ),
146 $ LDA )
147 CALL ZLACGV( N-K, A( K, K+1 ), LDA )
148 END IF
149 10 CONTINUE
150 ELSE
151 *
152 * Compute inv(L)*A*inv(L**H)
153 *
154 DO 20 K = 1, N
155 *
156 * Update the lower triangle of A(k:n,k:n)
157 *
158 AKK = A( K, K )
159 BKK = B( K, K )
160 AKK = AKK / BKK**2
161 A( K, K ) = AKK
162 IF( K.LT.N ) THEN
163 CALL ZDSCAL( N-K, ONE / BKK, A( K+1, K ), 1 )
164 CT = -HALF*AKK
165 CALL ZAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
166 CALL ZHER2( UPLO, N-K, -CONE, A( K+1, K ), 1,
167 $ B( K+1, K ), 1, A( K+1, K+1 ), LDA )
168 CALL ZAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
169 CALL ZTRSV( UPLO, 'No transpose', 'Non-unit', N-K,
170 $ B( K+1, K+1 ), LDB, A( K+1, K ), 1 )
171 END IF
172 20 CONTINUE
173 END IF
174 ELSE
175 IF( UPPER ) THEN
176 *
177 * Compute U*A*U**H
178 *
179 DO 30 K = 1, N
180 *
181 * Update the upper triangle of A(1:k,1:k)
182 *
183 AKK = A( K, K )
184 BKK = B( K, K )
185 CALL ZTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B,
186 $ LDB, A( 1, K ), 1 )
187 CT = HALF*AKK
188 CALL ZAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
189 CALL ZHER2( UPLO, K-1, CONE, A( 1, K ), 1, B( 1, K ), 1,
190 $ A, LDA )
191 CALL ZAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
192 CALL ZDSCAL( K-1, BKK, A( 1, K ), 1 )
193 A( K, K ) = AKK*BKK**2
194 30 CONTINUE
195 ELSE
196 *
197 * Compute L**H *A*L
198 *
199 DO 40 K = 1, N
200 *
201 * Update the lower triangle of A(1:k,1:k)
202 *
203 AKK = A( K, K )
204 BKK = B( K, K )
205 CALL ZLACGV( K-1, A( K, 1 ), LDA )
206 CALL ZTRMV( UPLO, 'Conjugate transpose', 'Non-unit', K-1,
207 $ B, LDB, A( K, 1 ), LDA )
208 CT = HALF*AKK
209 CALL ZLACGV( K-1, B( K, 1 ), LDB )
210 CALL ZAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
211 CALL ZHER2( UPLO, K-1, CONE, A( K, 1 ), LDA, B( K, 1 ),
212 $ LDB, A, LDA )
213 CALL ZAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
214 CALL ZLACGV( K-1, B( K, 1 ), LDB )
215 CALL ZDSCAL( K-1, BKK, A( K, 1 ), LDA )
216 CALL ZLACGV( K-1, A( K, 1 ), LDA )
217 A( K, K ) = AKK*BKK**2
218 40 CONTINUE
219 END IF
220 END IF
221 RETURN
222 *
223 * End of ZHEGS2
224 *
225 END