1 SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
2 $ LWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSYGV computes all the eigenvalues, and optionally, the eigenvectors
21 * of a real generalized symmetric-definite eigenproblem, of the form
22 * A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
23 * Here A and B are assumed to be symmetric and B is also
24 * positive definite.
25 *
26 * Arguments
27 * =========
28 *
29 * ITYPE (input) INTEGER
30 * Specifies the problem type to be solved:
31 * = 1: A*x = (lambda)*B*x
32 * = 2: A*B*x = (lambda)*x
33 * = 3: B*A*x = (lambda)*x
34 *
35 * JOBZ (input) CHARACTER*1
36 * = 'N': Compute eigenvalues only;
37 * = 'V': Compute eigenvalues and eigenvectors.
38 *
39 * UPLO (input) CHARACTER*1
40 * = 'U': Upper triangles of A and B are stored;
41 * = 'L': Lower triangles of A and B are stored.
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
48 * leading N-by-N upper triangular part of A contains the
49 * upper triangular part of the matrix A. If UPLO = 'L',
50 * the leading N-by-N lower triangular part of A contains
51 * the lower triangular part of the matrix A.
52 *
53 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
54 * matrix Z of eigenvectors. The eigenvectors are normalized
55 * as follows:
56 * if ITYPE = 1 or 2, Z**T*B*Z = I;
57 * if ITYPE = 3, Z**T*inv(B)*Z = I.
58 * If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
59 * or the lower triangle (if UPLO='L') of A, including the
60 * diagonal, is destroyed.
61 *
62 * LDA (input) INTEGER
63 * The leading dimension of the array A. LDA >= max(1,N).
64 *
65 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
66 * On entry, the symmetric positive definite matrix B.
67 * If UPLO = 'U', the leading N-by-N upper triangular part of B
68 * contains the upper triangular part of the matrix B.
69 * If UPLO = 'L', the leading N-by-N lower triangular part of B
70 * contains the lower triangular part of the matrix B.
71 *
72 * On exit, if INFO <= N, the part of B containing the matrix is
73 * overwritten by the triangular factor U or L from the Cholesky
74 * factorization B = U**T*U or B = L*L**T.
75 *
76 * LDB (input) INTEGER
77 * The leading dimension of the array B. LDB >= max(1,N).
78 *
79 * W (output) DOUBLE PRECISION array, dimension (N)
80 * If INFO = 0, the eigenvalues in ascending order.
81 *
82 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
83 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
84 *
85 * LWORK (input) INTEGER
86 * The length of the array WORK. LWORK >= max(1,3*N-1).
87 * For optimal efficiency, LWORK >= (NB+2)*N,
88 * where NB is the blocksize for DSYTRD returned by ILAENV.
89 *
90 * If LWORK = -1, then a workspace query is assumed; the routine
91 * only calculates the optimal size of the WORK array, returns
92 * this value as the first entry of the WORK array, and no error
93 * message related to LWORK is issued by XERBLA.
94 *
95 * INFO (output) INTEGER
96 * = 0: successful exit
97 * < 0: if INFO = -i, the i-th argument had an illegal value
98 * > 0: DPOTRF or DSYEV returned an error code:
99 * <= N: if INFO = i, DSYEV failed to converge;
100 * i off-diagonal elements of an intermediate
101 * tridiagonal form did not converge to zero;
102 * > N: if INFO = N + i, for 1 <= i <= N, then the leading
103 * minor of order i of B is not positive definite.
104 * The factorization of B could not be completed and
105 * no eigenvalues or eigenvectors were computed.
106 *
107 * =====================================================================
108 *
109 * .. Parameters ..
110 DOUBLE PRECISION ONE
111 PARAMETER ( ONE = 1.0D+0 )
112 * ..
113 * .. Local Scalars ..
114 LOGICAL LQUERY, UPPER, WANTZ
115 CHARACTER TRANS
116 INTEGER LWKMIN, LWKOPT, NB, NEIG
117 * ..
118 * .. External Functions ..
119 LOGICAL LSAME
120 INTEGER ILAENV
121 EXTERNAL LSAME, ILAENV
122 * ..
123 * .. External Subroutines ..
124 EXTERNAL DPOTRF, DSYEV, DSYGST, DTRMM, DTRSM, XERBLA
125 * ..
126 * .. Intrinsic Functions ..
127 INTRINSIC MAX
128 * ..
129 * .. Executable Statements ..
130 *
131 * Test the input parameters.
132 *
133 WANTZ = LSAME( JOBZ, 'V' )
134 UPPER = LSAME( UPLO, 'U' )
135 LQUERY = ( LWORK.EQ.-1 )
136 *
137 INFO = 0
138 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
139 INFO = -1
140 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
141 INFO = -2
142 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
143 INFO = -3
144 ELSE IF( N.LT.0 ) THEN
145 INFO = -4
146 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
147 INFO = -6
148 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
149 INFO = -8
150 END IF
151 *
152 IF( INFO.EQ.0 ) THEN
153 LWKMIN = MAX( 1, 3*N - 1 )
154 NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
155 LWKOPT = MAX( LWKMIN, ( NB + 2 )*N )
156 WORK( 1 ) = LWKOPT
157 *
158 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
159 INFO = -11
160 END IF
161 END IF
162 *
163 IF( INFO.NE.0 ) THEN
164 CALL XERBLA( 'DSYGV ', -INFO )
165 RETURN
166 ELSE IF( LQUERY ) THEN
167 RETURN
168 END IF
169 *
170 * Quick return if possible
171 *
172 IF( N.EQ.0 )
173 $ RETURN
174 *
175 * Form a Cholesky factorization of B.
176 *
177 CALL DPOTRF( UPLO, N, B, LDB, INFO )
178 IF( INFO.NE.0 ) THEN
179 INFO = N + INFO
180 RETURN
181 END IF
182 *
183 * Transform problem to standard eigenvalue problem and solve.
184 *
185 CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
186 CALL DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
187 *
188 IF( WANTZ ) THEN
189 *
190 * Backtransform eigenvectors to the original problem.
191 *
192 NEIG = N
193 IF( INFO.GT.0 )
194 $ NEIG = INFO - 1
195 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
196 *
197 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
198 * backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
199 *
200 IF( UPPER ) THEN
201 TRANS = 'N'
202 ELSE
203 TRANS = 'T'
204 END IF
205 *
206 CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
207 $ B, LDB, A, LDA )
208 *
209 ELSE IF( ITYPE.EQ.3 ) THEN
210 *
211 * For B*A*x=(lambda)*x;
212 * backtransform eigenvectors: x = L*y or U**T*y
213 *
214 IF( UPPER ) THEN
215 TRANS = 'T'
216 ELSE
217 TRANS = 'N'
218 END IF
219 *
220 CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
221 $ B, LDB, A, LDA )
222 END IF
223 END IF
224 *
225 WORK( 1 ) = LWKOPT
226 RETURN
227 *
228 * End of DSYGV
229 *
230 END
2 $ LWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSYGV computes all the eigenvalues, and optionally, the eigenvectors
21 * of a real generalized symmetric-definite eigenproblem, of the form
22 * A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
23 * Here A and B are assumed to be symmetric and B is also
24 * positive definite.
25 *
26 * Arguments
27 * =========
28 *
29 * ITYPE (input) INTEGER
30 * Specifies the problem type to be solved:
31 * = 1: A*x = (lambda)*B*x
32 * = 2: A*B*x = (lambda)*x
33 * = 3: B*A*x = (lambda)*x
34 *
35 * JOBZ (input) CHARACTER*1
36 * = 'N': Compute eigenvalues only;
37 * = 'V': Compute eigenvalues and eigenvectors.
38 *
39 * UPLO (input) CHARACTER*1
40 * = 'U': Upper triangles of A and B are stored;
41 * = 'L': Lower triangles of A and B are stored.
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
48 * leading N-by-N upper triangular part of A contains the
49 * upper triangular part of the matrix A. If UPLO = 'L',
50 * the leading N-by-N lower triangular part of A contains
51 * the lower triangular part of the matrix A.
52 *
53 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
54 * matrix Z of eigenvectors. The eigenvectors are normalized
55 * as follows:
56 * if ITYPE = 1 or 2, Z**T*B*Z = I;
57 * if ITYPE = 3, Z**T*inv(B)*Z = I.
58 * If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
59 * or the lower triangle (if UPLO='L') of A, including the
60 * diagonal, is destroyed.
61 *
62 * LDA (input) INTEGER
63 * The leading dimension of the array A. LDA >= max(1,N).
64 *
65 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
66 * On entry, the symmetric positive definite matrix B.
67 * If UPLO = 'U', the leading N-by-N upper triangular part of B
68 * contains the upper triangular part of the matrix B.
69 * If UPLO = 'L', the leading N-by-N lower triangular part of B
70 * contains the lower triangular part of the matrix B.
71 *
72 * On exit, if INFO <= N, the part of B containing the matrix is
73 * overwritten by the triangular factor U or L from the Cholesky
74 * factorization B = U**T*U or B = L*L**T.
75 *
76 * LDB (input) INTEGER
77 * The leading dimension of the array B. LDB >= max(1,N).
78 *
79 * W (output) DOUBLE PRECISION array, dimension (N)
80 * If INFO = 0, the eigenvalues in ascending order.
81 *
82 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
83 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
84 *
85 * LWORK (input) INTEGER
86 * The length of the array WORK. LWORK >= max(1,3*N-1).
87 * For optimal efficiency, LWORK >= (NB+2)*N,
88 * where NB is the blocksize for DSYTRD returned by ILAENV.
89 *
90 * If LWORK = -1, then a workspace query is assumed; the routine
91 * only calculates the optimal size of the WORK array, returns
92 * this value as the first entry of the WORK array, and no error
93 * message related to LWORK is issued by XERBLA.
94 *
95 * INFO (output) INTEGER
96 * = 0: successful exit
97 * < 0: if INFO = -i, the i-th argument had an illegal value
98 * > 0: DPOTRF or DSYEV returned an error code:
99 * <= N: if INFO = i, DSYEV failed to converge;
100 * i off-diagonal elements of an intermediate
101 * tridiagonal form did not converge to zero;
102 * > N: if INFO = N + i, for 1 <= i <= N, then the leading
103 * minor of order i of B is not positive definite.
104 * The factorization of B could not be completed and
105 * no eigenvalues or eigenvectors were computed.
106 *
107 * =====================================================================
108 *
109 * .. Parameters ..
110 DOUBLE PRECISION ONE
111 PARAMETER ( ONE = 1.0D+0 )
112 * ..
113 * .. Local Scalars ..
114 LOGICAL LQUERY, UPPER, WANTZ
115 CHARACTER TRANS
116 INTEGER LWKMIN, LWKOPT, NB, NEIG
117 * ..
118 * .. External Functions ..
119 LOGICAL LSAME
120 INTEGER ILAENV
121 EXTERNAL LSAME, ILAENV
122 * ..
123 * .. External Subroutines ..
124 EXTERNAL DPOTRF, DSYEV, DSYGST, DTRMM, DTRSM, XERBLA
125 * ..
126 * .. Intrinsic Functions ..
127 INTRINSIC MAX
128 * ..
129 * .. Executable Statements ..
130 *
131 * Test the input parameters.
132 *
133 WANTZ = LSAME( JOBZ, 'V' )
134 UPPER = LSAME( UPLO, 'U' )
135 LQUERY = ( LWORK.EQ.-1 )
136 *
137 INFO = 0
138 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
139 INFO = -1
140 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
141 INFO = -2
142 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
143 INFO = -3
144 ELSE IF( N.LT.0 ) THEN
145 INFO = -4
146 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
147 INFO = -6
148 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
149 INFO = -8
150 END IF
151 *
152 IF( INFO.EQ.0 ) THEN
153 LWKMIN = MAX( 1, 3*N - 1 )
154 NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
155 LWKOPT = MAX( LWKMIN, ( NB + 2 )*N )
156 WORK( 1 ) = LWKOPT
157 *
158 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
159 INFO = -11
160 END IF
161 END IF
162 *
163 IF( INFO.NE.0 ) THEN
164 CALL XERBLA( 'DSYGV ', -INFO )
165 RETURN
166 ELSE IF( LQUERY ) THEN
167 RETURN
168 END IF
169 *
170 * Quick return if possible
171 *
172 IF( N.EQ.0 )
173 $ RETURN
174 *
175 * Form a Cholesky factorization of B.
176 *
177 CALL DPOTRF( UPLO, N, B, LDB, INFO )
178 IF( INFO.NE.0 ) THEN
179 INFO = N + INFO
180 RETURN
181 END IF
182 *
183 * Transform problem to standard eigenvalue problem and solve.
184 *
185 CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
186 CALL DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
187 *
188 IF( WANTZ ) THEN
189 *
190 * Backtransform eigenvectors to the original problem.
191 *
192 NEIG = N
193 IF( INFO.GT.0 )
194 $ NEIG = INFO - 1
195 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
196 *
197 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
198 * backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
199 *
200 IF( UPPER ) THEN
201 TRANS = 'N'
202 ELSE
203 TRANS = 'T'
204 END IF
205 *
206 CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
207 $ B, LDB, A, LDA )
208 *
209 ELSE IF( ITYPE.EQ.3 ) THEN
210 *
211 * For B*A*x=(lambda)*x;
212 * backtransform eigenvectors: x = L*y or U**T*y
213 *
214 IF( UPPER ) THEN
215 TRANS = 'T'
216 ELSE
217 TRANS = 'N'
218 END IF
219 *
220 CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
221 $ B, LDB, A, LDA )
222 END IF
223 END IF
224 *
225 WORK( 1 ) = LWKOPT
226 RETURN
227 *
228 * End of DSYGV
229 *
230 END