1 SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
2 *
3 * -- LAPACK driver routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER JOBZ, UPLO
10 INTEGER INFO, LDA, LWORK, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSYEV computes all eigenvalues and, optionally, eigenvectors of a
20 * real symmetric matrix A.
21 *
22 * Arguments
23 * =========
24 *
25 * JOBZ (input) CHARACTER*1
26 * = 'N': Compute eigenvalues only;
27 * = 'V': Compute eigenvalues and eigenvectors.
28 *
29 * UPLO (input) CHARACTER*1
30 * = 'U': Upper triangle of A is stored;
31 * = 'L': Lower triangle of A is stored.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
37 * On entry, the symmetric matrix A. If UPLO = 'U', the
38 * leading N-by-N upper triangular part of A contains the
39 * upper triangular part of the matrix A. If UPLO = 'L',
40 * the leading N-by-N lower triangular part of A contains
41 * the lower triangular part of the matrix A.
42 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
43 * orthonormal eigenvectors of the matrix A.
44 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
45 * or the upper triangle (if UPLO='U') of A, including the
46 * diagonal, is destroyed.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,N).
50 *
51 * W (output) DOUBLE PRECISION array, dimension (N)
52 * If INFO = 0, the eigenvalues in ascending order.
53 *
54 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
55 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
56 *
57 * LWORK (input) INTEGER
58 * The length of the array WORK. LWORK >= max(1,3*N-1).
59 * For optimal efficiency, LWORK >= (NB+2)*N,
60 * where NB is the blocksize for DSYTRD returned by ILAENV.
61 *
62 * If LWORK = -1, then a workspace query is assumed; the routine
63 * only calculates the optimal size of the WORK array, returns
64 * this value as the first entry of the WORK array, and no error
65 * message related to LWORK is issued by XERBLA.
66 *
67 * INFO (output) INTEGER
68 * = 0: successful exit
69 * < 0: if INFO = -i, the i-th argument had an illegal value
70 * > 0: if INFO = i, the algorithm failed to converge; i
71 * off-diagonal elements of an intermediate tridiagonal
72 * form did not converge to zero.
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77 DOUBLE PRECISION ZERO, ONE
78 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
79 * ..
80 * .. Local Scalars ..
81 LOGICAL LOWER, LQUERY, WANTZ
82 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
83 $ LLWORK, LWKOPT, NB
84 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
85 $ SMLNUM
86 * ..
87 * .. External Functions ..
88 LOGICAL LSAME
89 INTEGER ILAENV
90 DOUBLE PRECISION DLAMCH, DLANSY
91 EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD,
95 $ XERBLA
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC MAX, SQRT
99 * ..
100 * .. Executable Statements ..
101 *
102 * Test the input parameters.
103 *
104 WANTZ = LSAME( JOBZ, 'V' )
105 LOWER = LSAME( UPLO, 'L' )
106 LQUERY = ( LWORK.EQ.-1 )
107 *
108 INFO = 0
109 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
110 INFO = -1
111 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
112 INFO = -2
113 ELSE IF( N.LT.0 ) THEN
114 INFO = -3
115 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
116 INFO = -5
117 END IF
118 *
119 IF( INFO.EQ.0 ) THEN
120 NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
121 LWKOPT = MAX( 1, ( NB+2 )*N )
122 WORK( 1 ) = LWKOPT
123 *
124 IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY )
125 $ INFO = -8
126 END IF
127 *
128 IF( INFO.NE.0 ) THEN
129 CALL XERBLA( 'DSYEV ', -INFO )
130 RETURN
131 ELSE IF( LQUERY ) THEN
132 RETURN
133 END IF
134 *
135 * Quick return if possible
136 *
137 IF( N.EQ.0 ) THEN
138 RETURN
139 END IF
140 *
141 IF( N.EQ.1 ) THEN
142 W( 1 ) = A( 1, 1 )
143 WORK( 1 ) = 2
144 IF( WANTZ )
145 $ A( 1, 1 ) = ONE
146 RETURN
147 END IF
148 *
149 * Get machine constants.
150 *
151 SAFMIN = DLAMCH( 'Safe minimum' )
152 EPS = DLAMCH( 'Precision' )
153 SMLNUM = SAFMIN / EPS
154 BIGNUM = ONE / SMLNUM
155 RMIN = SQRT( SMLNUM )
156 RMAX = SQRT( BIGNUM )
157 *
158 * Scale matrix to allowable range, if necessary.
159 *
160 ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
161 ISCALE = 0
162 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
163 ISCALE = 1
164 SIGMA = RMIN / ANRM
165 ELSE IF( ANRM.GT.RMAX ) THEN
166 ISCALE = 1
167 SIGMA = RMAX / ANRM
168 END IF
169 IF( ISCALE.EQ.1 )
170 $ CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
171 *
172 * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
173 *
174 INDE = 1
175 INDTAU = INDE + N
176 INDWRK = INDTAU + N
177 LLWORK = LWORK - INDWRK + 1
178 CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
179 $ WORK( INDWRK ), LLWORK, IINFO )
180 *
181 * For eigenvalues only, call DSTERF. For eigenvectors, first call
182 * DORGTR to generate the orthogonal matrix, then call DSTEQR.
183 *
184 IF( .NOT.WANTZ ) THEN
185 CALL DSTERF( N, W, WORK( INDE ), INFO )
186 ELSE
187 CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
188 $ LLWORK, IINFO )
189 CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ),
190 $ INFO )
191 END IF
192 *
193 * If matrix was scaled, then rescale eigenvalues appropriately.
194 *
195 IF( ISCALE.EQ.1 ) THEN
196 IF( INFO.EQ.0 ) THEN
197 IMAX = N
198 ELSE
199 IMAX = INFO - 1
200 END IF
201 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
202 END IF
203 *
204 * Set WORK(1) to optimal workspace size.
205 *
206 WORK( 1 ) = LWKOPT
207 *
208 RETURN
209 *
210 * End of DSYEV
211 *
212 END
2 *
3 * -- LAPACK driver routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER JOBZ, UPLO
10 INTEGER INFO, LDA, LWORK, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSYEV computes all eigenvalues and, optionally, eigenvectors of a
20 * real symmetric matrix A.
21 *
22 * Arguments
23 * =========
24 *
25 * JOBZ (input) CHARACTER*1
26 * = 'N': Compute eigenvalues only;
27 * = 'V': Compute eigenvalues and eigenvectors.
28 *
29 * UPLO (input) CHARACTER*1
30 * = 'U': Upper triangle of A is stored;
31 * = 'L': Lower triangle of A is stored.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
37 * On entry, the symmetric matrix A. If UPLO = 'U', the
38 * leading N-by-N upper triangular part of A contains the
39 * upper triangular part of the matrix A. If UPLO = 'L',
40 * the leading N-by-N lower triangular part of A contains
41 * the lower triangular part of the matrix A.
42 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
43 * orthonormal eigenvectors of the matrix A.
44 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
45 * or the upper triangle (if UPLO='U') of A, including the
46 * diagonal, is destroyed.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,N).
50 *
51 * W (output) DOUBLE PRECISION array, dimension (N)
52 * If INFO = 0, the eigenvalues in ascending order.
53 *
54 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
55 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
56 *
57 * LWORK (input) INTEGER
58 * The length of the array WORK. LWORK >= max(1,3*N-1).
59 * For optimal efficiency, LWORK >= (NB+2)*N,
60 * where NB is the blocksize for DSYTRD returned by ILAENV.
61 *
62 * If LWORK = -1, then a workspace query is assumed; the routine
63 * only calculates the optimal size of the WORK array, returns
64 * this value as the first entry of the WORK array, and no error
65 * message related to LWORK is issued by XERBLA.
66 *
67 * INFO (output) INTEGER
68 * = 0: successful exit
69 * < 0: if INFO = -i, the i-th argument had an illegal value
70 * > 0: if INFO = i, the algorithm failed to converge; i
71 * off-diagonal elements of an intermediate tridiagonal
72 * form did not converge to zero.
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77 DOUBLE PRECISION ZERO, ONE
78 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
79 * ..
80 * .. Local Scalars ..
81 LOGICAL LOWER, LQUERY, WANTZ
82 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
83 $ LLWORK, LWKOPT, NB
84 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
85 $ SMLNUM
86 * ..
87 * .. External Functions ..
88 LOGICAL LSAME
89 INTEGER ILAENV
90 DOUBLE PRECISION DLAMCH, DLANSY
91 EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD,
95 $ XERBLA
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC MAX, SQRT
99 * ..
100 * .. Executable Statements ..
101 *
102 * Test the input parameters.
103 *
104 WANTZ = LSAME( JOBZ, 'V' )
105 LOWER = LSAME( UPLO, 'L' )
106 LQUERY = ( LWORK.EQ.-1 )
107 *
108 INFO = 0
109 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
110 INFO = -1
111 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
112 INFO = -2
113 ELSE IF( N.LT.0 ) THEN
114 INFO = -3
115 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
116 INFO = -5
117 END IF
118 *
119 IF( INFO.EQ.0 ) THEN
120 NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
121 LWKOPT = MAX( 1, ( NB+2 )*N )
122 WORK( 1 ) = LWKOPT
123 *
124 IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY )
125 $ INFO = -8
126 END IF
127 *
128 IF( INFO.NE.0 ) THEN
129 CALL XERBLA( 'DSYEV ', -INFO )
130 RETURN
131 ELSE IF( LQUERY ) THEN
132 RETURN
133 END IF
134 *
135 * Quick return if possible
136 *
137 IF( N.EQ.0 ) THEN
138 RETURN
139 END IF
140 *
141 IF( N.EQ.1 ) THEN
142 W( 1 ) = A( 1, 1 )
143 WORK( 1 ) = 2
144 IF( WANTZ )
145 $ A( 1, 1 ) = ONE
146 RETURN
147 END IF
148 *
149 * Get machine constants.
150 *
151 SAFMIN = DLAMCH( 'Safe minimum' )
152 EPS = DLAMCH( 'Precision' )
153 SMLNUM = SAFMIN / EPS
154 BIGNUM = ONE / SMLNUM
155 RMIN = SQRT( SMLNUM )
156 RMAX = SQRT( BIGNUM )
157 *
158 * Scale matrix to allowable range, if necessary.
159 *
160 ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
161 ISCALE = 0
162 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
163 ISCALE = 1
164 SIGMA = RMIN / ANRM
165 ELSE IF( ANRM.GT.RMAX ) THEN
166 ISCALE = 1
167 SIGMA = RMAX / ANRM
168 END IF
169 IF( ISCALE.EQ.1 )
170 $ CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
171 *
172 * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
173 *
174 INDE = 1
175 INDTAU = INDE + N
176 INDWRK = INDTAU + N
177 LLWORK = LWORK - INDWRK + 1
178 CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
179 $ WORK( INDWRK ), LLWORK, IINFO )
180 *
181 * For eigenvalues only, call DSTERF. For eigenvectors, first call
182 * DORGTR to generate the orthogonal matrix, then call DSTEQR.
183 *
184 IF( .NOT.WANTZ ) THEN
185 CALL DSTERF( N, W, WORK( INDE ), INFO )
186 ELSE
187 CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
188 $ LLWORK, IINFO )
189 CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ),
190 $ INFO )
191 END IF
192 *
193 * If matrix was scaled, then rescale eigenvalues appropriately.
194 *
195 IF( ISCALE.EQ.1 ) THEN
196 IF( INFO.EQ.0 ) THEN
197 IMAX = N
198 ELSE
199 IMAX = INFO - 1
200 END IF
201 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
202 END IF
203 *
204 * Set WORK(1) to optimal workspace size.
205 *
206 WORK( 1 ) = LWKOPT
207 *
208 RETURN
209 *
210 * End of DSYEV
211 *
212 END