1 SUBROUTINE DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK,
2 $ LIWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, LDA, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSYEVD computes all eigenvalues and, optionally, eigenvectors of a
22 * real symmetric matrix A. If eigenvectors are desired, it uses a
23 * divide and conquer algorithm.
24 *
25 * The divide and conquer algorithm makes very mild assumptions about
26 * floating point arithmetic. It will work on machines with a guard
27 * digit in add/subtract, or on those binary machines without guard
28 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
29 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
30 * without guard digits, but we know of none.
31 *
32 * Because of large use of BLAS of level 3, DSYEVD needs N**2 more
33 * workspace than DSYEVX.
34 *
35 * Arguments
36 * =========
37 *
38 * JOBZ (input) CHARACTER*1
39 * = 'N': Compute eigenvalues only;
40 * = 'V': Compute eigenvalues and eigenvectors.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangle of A is stored;
44 * = 'L': Lower triangle of A is stored.
45 *
46 * N (input) INTEGER
47 * The order of the matrix A. N >= 0.
48 *
49 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
50 * On entry, the symmetric matrix A. If UPLO = 'U', the
51 * leading N-by-N upper triangular part of A contains the
52 * upper triangular part of the matrix A. If UPLO = 'L',
53 * the leading N-by-N lower triangular part of A contains
54 * the lower triangular part of the matrix A.
55 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
56 * orthonormal eigenvectors of the matrix A.
57 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
58 * or the upper triangle (if UPLO='U') of A, including the
59 * diagonal, is destroyed.
60 *
61 * LDA (input) INTEGER
62 * The leading dimension of the array A. LDA >= max(1,N).
63 *
64 * W (output) DOUBLE PRECISION array, dimension (N)
65 * If INFO = 0, the eigenvalues in ascending order.
66 *
67 * WORK (workspace/output) DOUBLE PRECISION array,
68 * dimension (LWORK)
69 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
70 *
71 * LWORK (input) INTEGER
72 * The dimension of the array WORK.
73 * If N <= 1, LWORK must be at least 1.
74 * If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
75 * If JOBZ = 'V' and N > 1, LWORK must be at least
76 * 1 + 6*N + 2*N**2.
77 *
78 * If LWORK = -1, then a workspace query is assumed; the routine
79 * only calculates the optimal sizes of the WORK and IWORK
80 * arrays, returns these values as the first entries of the WORK
81 * and IWORK arrays, and no error message related to LWORK or
82 * LIWORK is issued by XERBLA.
83 *
84 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
85 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
86 *
87 * LIWORK (input) INTEGER
88 * The dimension of the array IWORK.
89 * If N <= 1, LIWORK must be at least 1.
90 * If JOBZ = 'N' and N > 1, LIWORK must be at least 1.
91 * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
92 *
93 * If LIWORK = -1, then a workspace query is assumed; the
94 * routine only calculates the optimal sizes of the WORK and
95 * IWORK arrays, returns these values as the first entries of
96 * the WORK and IWORK arrays, and no error message related to
97 * LWORK or LIWORK is issued by XERBLA.
98 *
99 * INFO (output) INTEGER
100 * = 0: successful exit
101 * < 0: if INFO = -i, the i-th argument had an illegal value
102 * > 0: if INFO = i and JOBZ = 'N', then the algorithm failed
103 * to converge; i off-diagonal elements of an intermediate
104 * tridiagonal form did not converge to zero;
105 * if INFO = i and JOBZ = 'V', then the algorithm failed
106 * to compute an eigenvalue while working on the submatrix
107 * lying in rows and columns INFO/(N+1) through
108 * mod(INFO,N+1).
109 *
110 * Further Details
111 * ===============
112 *
113 * Based on contributions by
114 * Jeff Rutter, Computer Science Division, University of California
115 * at Berkeley, USA
116 * Modified by Francoise Tisseur, University of Tennessee.
117 *
118 * Modified description of INFO. Sven, 16 Feb 05.
119 * =====================================================================
120 *
121 * .. Parameters ..
122 DOUBLE PRECISION ZERO, ONE
123 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
124 * ..
125 * .. Local Scalars ..
126 *
127 LOGICAL LOWER, LQUERY, WANTZ
128 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
129 $ LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN
130 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
131 $ SMLNUM
132 * ..
133 * .. External Functions ..
134 LOGICAL LSAME
135 INTEGER ILAENV
136 DOUBLE PRECISION DLAMCH, DLANSY
137 EXTERNAL LSAME, DLAMCH, DLANSY, ILAENV
138 * ..
139 * .. External Subroutines ..
140 EXTERNAL DLACPY, DLASCL, DORMTR, DSCAL, DSTEDC, DSTERF,
141 $ DSYTRD, XERBLA
142 * ..
143 * .. Intrinsic Functions ..
144 INTRINSIC MAX, SQRT
145 * ..
146 * .. Executable Statements ..
147 *
148 * Test the input parameters.
149 *
150 WANTZ = LSAME( JOBZ, 'V' )
151 LOWER = LSAME( UPLO, 'L' )
152 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
153 *
154 INFO = 0
155 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
156 INFO = -1
157 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
158 INFO = -2
159 ELSE IF( N.LT.0 ) THEN
160 INFO = -3
161 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
162 INFO = -5
163 END IF
164 *
165 IF( INFO.EQ.0 ) THEN
166 IF( N.LE.1 ) THEN
167 LIWMIN = 1
168 LWMIN = 1
169 LOPT = LWMIN
170 LIOPT = LIWMIN
171 ELSE
172 IF( WANTZ ) THEN
173 LIWMIN = 3 + 5*N
174 LWMIN = 1 + 6*N + 2*N**2
175 ELSE
176 LIWMIN = 1
177 LWMIN = 2*N + 1
178 END IF
179 LOPT = MAX( LWMIN, 2*N +
180 $ ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
181 LIOPT = LIWMIN
182 END IF
183 WORK( 1 ) = LOPT
184 IWORK( 1 ) = LIOPT
185 *
186 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
187 INFO = -8
188 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
189 INFO = -10
190 END IF
191 END IF
192 *
193 IF( INFO.NE.0 ) THEN
194 CALL XERBLA( 'DSYEVD', -INFO )
195 RETURN
196 ELSE IF( LQUERY ) THEN
197 RETURN
198 END IF
199 *
200 * Quick return if possible
201 *
202 IF( N.EQ.0 )
203 $ RETURN
204 *
205 IF( N.EQ.1 ) THEN
206 W( 1 ) = A( 1, 1 )
207 IF( WANTZ )
208 $ A( 1, 1 ) = ONE
209 RETURN
210 END IF
211 *
212 * Get machine constants.
213 *
214 SAFMIN = DLAMCH( 'Safe minimum' )
215 EPS = DLAMCH( 'Precision' )
216 SMLNUM = SAFMIN / EPS
217 BIGNUM = ONE / SMLNUM
218 RMIN = SQRT( SMLNUM )
219 RMAX = SQRT( BIGNUM )
220 *
221 * Scale matrix to allowable range, if necessary.
222 *
223 ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
224 ISCALE = 0
225 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
226 ISCALE = 1
227 SIGMA = RMIN / ANRM
228 ELSE IF( ANRM.GT.RMAX ) THEN
229 ISCALE = 1
230 SIGMA = RMAX / ANRM
231 END IF
232 IF( ISCALE.EQ.1 )
233 $ CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
234 *
235 * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
236 *
237 INDE = 1
238 INDTAU = INDE + N
239 INDWRK = INDTAU + N
240 LLWORK = LWORK - INDWRK + 1
241 INDWK2 = INDWRK + N*N
242 LLWRK2 = LWORK - INDWK2 + 1
243 *
244 CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
245 $ WORK( INDWRK ), LLWORK, IINFO )
246 LOPT = 2*N + WORK( INDWRK )
247 *
248 * For eigenvalues only, call DSTERF. For eigenvectors, first call
249 * DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
250 * tridiagonal matrix, then call DORMTR to multiply it by the
251 * Householder transformations stored in A.
252 *
253 IF( .NOT.WANTZ ) THEN
254 CALL DSTERF( N, W, WORK( INDE ), INFO )
255 ELSE
256 CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
257 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
258 CALL DORMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
259 $ WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
260 CALL DLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
261 LOPT = MAX( LOPT, 1+6*N+2*N**2 )
262 END IF
263 *
264 * If matrix was scaled, then rescale eigenvalues appropriately.
265 *
266 IF( ISCALE.EQ.1 )
267 $ CALL DSCAL( N, ONE / SIGMA, W, 1 )
268 *
269 WORK( 1 ) = LOPT
270 IWORK( 1 ) = LIOPT
271 *
272 RETURN
273 *
274 * End of DSYEVD
275 *
276 END
2 $ LIWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, LDA, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSYEVD computes all eigenvalues and, optionally, eigenvectors of a
22 * real symmetric matrix A. If eigenvectors are desired, it uses a
23 * divide and conquer algorithm.
24 *
25 * The divide and conquer algorithm makes very mild assumptions about
26 * floating point arithmetic. It will work on machines with a guard
27 * digit in add/subtract, or on those binary machines without guard
28 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
29 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
30 * without guard digits, but we know of none.
31 *
32 * Because of large use of BLAS of level 3, DSYEVD needs N**2 more
33 * workspace than DSYEVX.
34 *
35 * Arguments
36 * =========
37 *
38 * JOBZ (input) CHARACTER*1
39 * = 'N': Compute eigenvalues only;
40 * = 'V': Compute eigenvalues and eigenvectors.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangle of A is stored;
44 * = 'L': Lower triangle of A is stored.
45 *
46 * N (input) INTEGER
47 * The order of the matrix A. N >= 0.
48 *
49 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
50 * On entry, the symmetric matrix A. If UPLO = 'U', the
51 * leading N-by-N upper triangular part of A contains the
52 * upper triangular part of the matrix A. If UPLO = 'L',
53 * the leading N-by-N lower triangular part of A contains
54 * the lower triangular part of the matrix A.
55 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
56 * orthonormal eigenvectors of the matrix A.
57 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
58 * or the upper triangle (if UPLO='U') of A, including the
59 * diagonal, is destroyed.
60 *
61 * LDA (input) INTEGER
62 * The leading dimension of the array A. LDA >= max(1,N).
63 *
64 * W (output) DOUBLE PRECISION array, dimension (N)
65 * If INFO = 0, the eigenvalues in ascending order.
66 *
67 * WORK (workspace/output) DOUBLE PRECISION array,
68 * dimension (LWORK)
69 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
70 *
71 * LWORK (input) INTEGER
72 * The dimension of the array WORK.
73 * If N <= 1, LWORK must be at least 1.
74 * If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
75 * If JOBZ = 'V' and N > 1, LWORK must be at least
76 * 1 + 6*N + 2*N**2.
77 *
78 * If LWORK = -1, then a workspace query is assumed; the routine
79 * only calculates the optimal sizes of the WORK and IWORK
80 * arrays, returns these values as the first entries of the WORK
81 * and IWORK arrays, and no error message related to LWORK or
82 * LIWORK is issued by XERBLA.
83 *
84 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
85 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
86 *
87 * LIWORK (input) INTEGER
88 * The dimension of the array IWORK.
89 * If N <= 1, LIWORK must be at least 1.
90 * If JOBZ = 'N' and N > 1, LIWORK must be at least 1.
91 * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
92 *
93 * If LIWORK = -1, then a workspace query is assumed; the
94 * routine only calculates the optimal sizes of the WORK and
95 * IWORK arrays, returns these values as the first entries of
96 * the WORK and IWORK arrays, and no error message related to
97 * LWORK or LIWORK is issued by XERBLA.
98 *
99 * INFO (output) INTEGER
100 * = 0: successful exit
101 * < 0: if INFO = -i, the i-th argument had an illegal value
102 * > 0: if INFO = i and JOBZ = 'N', then the algorithm failed
103 * to converge; i off-diagonal elements of an intermediate
104 * tridiagonal form did not converge to zero;
105 * if INFO = i and JOBZ = 'V', then the algorithm failed
106 * to compute an eigenvalue while working on the submatrix
107 * lying in rows and columns INFO/(N+1) through
108 * mod(INFO,N+1).
109 *
110 * Further Details
111 * ===============
112 *
113 * Based on contributions by
114 * Jeff Rutter, Computer Science Division, University of California
115 * at Berkeley, USA
116 * Modified by Francoise Tisseur, University of Tennessee.
117 *
118 * Modified description of INFO. Sven, 16 Feb 05.
119 * =====================================================================
120 *
121 * .. Parameters ..
122 DOUBLE PRECISION ZERO, ONE
123 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
124 * ..
125 * .. Local Scalars ..
126 *
127 LOGICAL LOWER, LQUERY, WANTZ
128 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
129 $ LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN
130 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
131 $ SMLNUM
132 * ..
133 * .. External Functions ..
134 LOGICAL LSAME
135 INTEGER ILAENV
136 DOUBLE PRECISION DLAMCH, DLANSY
137 EXTERNAL LSAME, DLAMCH, DLANSY, ILAENV
138 * ..
139 * .. External Subroutines ..
140 EXTERNAL DLACPY, DLASCL, DORMTR, DSCAL, DSTEDC, DSTERF,
141 $ DSYTRD, XERBLA
142 * ..
143 * .. Intrinsic Functions ..
144 INTRINSIC MAX, SQRT
145 * ..
146 * .. Executable Statements ..
147 *
148 * Test the input parameters.
149 *
150 WANTZ = LSAME( JOBZ, 'V' )
151 LOWER = LSAME( UPLO, 'L' )
152 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
153 *
154 INFO = 0
155 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
156 INFO = -1
157 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
158 INFO = -2
159 ELSE IF( N.LT.0 ) THEN
160 INFO = -3
161 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
162 INFO = -5
163 END IF
164 *
165 IF( INFO.EQ.0 ) THEN
166 IF( N.LE.1 ) THEN
167 LIWMIN = 1
168 LWMIN = 1
169 LOPT = LWMIN
170 LIOPT = LIWMIN
171 ELSE
172 IF( WANTZ ) THEN
173 LIWMIN = 3 + 5*N
174 LWMIN = 1 + 6*N + 2*N**2
175 ELSE
176 LIWMIN = 1
177 LWMIN = 2*N + 1
178 END IF
179 LOPT = MAX( LWMIN, 2*N +
180 $ ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
181 LIOPT = LIWMIN
182 END IF
183 WORK( 1 ) = LOPT
184 IWORK( 1 ) = LIOPT
185 *
186 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
187 INFO = -8
188 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
189 INFO = -10
190 END IF
191 END IF
192 *
193 IF( INFO.NE.0 ) THEN
194 CALL XERBLA( 'DSYEVD', -INFO )
195 RETURN
196 ELSE IF( LQUERY ) THEN
197 RETURN
198 END IF
199 *
200 * Quick return if possible
201 *
202 IF( N.EQ.0 )
203 $ RETURN
204 *
205 IF( N.EQ.1 ) THEN
206 W( 1 ) = A( 1, 1 )
207 IF( WANTZ )
208 $ A( 1, 1 ) = ONE
209 RETURN
210 END IF
211 *
212 * Get machine constants.
213 *
214 SAFMIN = DLAMCH( 'Safe minimum' )
215 EPS = DLAMCH( 'Precision' )
216 SMLNUM = SAFMIN / EPS
217 BIGNUM = ONE / SMLNUM
218 RMIN = SQRT( SMLNUM )
219 RMAX = SQRT( BIGNUM )
220 *
221 * Scale matrix to allowable range, if necessary.
222 *
223 ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
224 ISCALE = 0
225 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
226 ISCALE = 1
227 SIGMA = RMIN / ANRM
228 ELSE IF( ANRM.GT.RMAX ) THEN
229 ISCALE = 1
230 SIGMA = RMAX / ANRM
231 END IF
232 IF( ISCALE.EQ.1 )
233 $ CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
234 *
235 * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
236 *
237 INDE = 1
238 INDTAU = INDE + N
239 INDWRK = INDTAU + N
240 LLWORK = LWORK - INDWRK + 1
241 INDWK2 = INDWRK + N*N
242 LLWRK2 = LWORK - INDWK2 + 1
243 *
244 CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
245 $ WORK( INDWRK ), LLWORK, IINFO )
246 LOPT = 2*N + WORK( INDWRK )
247 *
248 * For eigenvalues only, call DSTERF. For eigenvectors, first call
249 * DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
250 * tridiagonal matrix, then call DORMTR to multiply it by the
251 * Householder transformations stored in A.
252 *
253 IF( .NOT.WANTZ ) THEN
254 CALL DSTERF( N, W, WORK( INDE ), INFO )
255 ELSE
256 CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
257 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
258 CALL DORMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
259 $ WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
260 CALL DLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
261 LOPT = MAX( LOPT, 1+6*N+2*N**2 )
262 END IF
263 *
264 * If matrix was scaled, then rescale eigenvalues appropriately.
265 *
266 IF( ISCALE.EQ.1 )
267 $ CALL DSCAL( N, ONE / SIGMA, W, 1 )
268 *
269 WORK( 1 ) = LOPT
270 IWORK( 1 ) = LIOPT
271 *
272 RETURN
273 *
274 * End of DSYEVD
275 *
276 END