| blob | 9d316112f8f741e1529a6a3c32eb1e6b8febe70e |
1 SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
2 !
3 ! -- LAPACK driver routine (version 3.1) --
4 ! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 ! November 2006
6 !
7 ! .. Scalar Arguments ..
8 CHARACTER JOBZ, UPLO
9 INTEGER INFO, LDA, LWORK, N
10 ! ..
11 ! .. Array Arguments ..
12 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
13 ! ..
14 !
15 ! Purpose
16 ! =======
17 !
18 ! DSYEV computes all eigenvalues and, optionally, eigenvectors of a
19 ! real symmetric matrix A.
20 !
21 ! Arguments
22 ! =========
23 !
24 ! JOBZ (input) CHARACTER*1
25 ! = 'N': Compute eigenvalues only;
26 ! = 'V': Compute eigenvalues and eigenvectors.
27 !
28 ! UPLO (input) CHARACTER*1
29 ! = 'U': Upper triangle of A is stored;
30 ! = 'L': Lower triangle of A is stored.
31 !
32 ! N (input) INTEGER
33 ! The order of the matrix A. N >= 0.
34 !
35 ! A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
36 ! On entry, the symmetric matrix A. If UPLO = 'U', the
37 ! leading N-by-N upper triangular part of A contains the
38 ! upper triangular part of the matrix A. If UPLO = 'L',
39 ! the leading N-by-N lower triangular part of A contains
40 ! the lower triangular part of the matrix A.
41 ! On exit, if JOBZ = 'V', then if INFO = 0, A contains the
42 ! orthonormal eigenvectors of the matrix A.
43 ! If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
44 ! or the upper triangle (if UPLO='U') of A, including the
45 ! diagonal, is destroyed.
46 !
47 ! LDA (input) INTEGER
48 ! The leading dimension of the array A. LDA >= max(1,N).
49 !
50 ! W (output) DOUBLE PRECISION array, dimension (N)
51 ! If INFO = 0, the eigenvalues in ascending order.
52 !
53 ! WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
54 ! On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
55 !
56 ! LWORK (input) INTEGER
57 ! The length of the array WORK. LWORK >= max(1,3*N-1).
58 ! For optimal efficiency, LWORK >= (NB+2)*N,
59 ! where NB is the blocksize for DSYTRD returned by ILAENV.
60 !
61 ! If LWORK = -1, then a workspace query is assumed; the routine
62 ! only calculates the optimal size of the WORK array, returns
63 ! this value as the first entry of the WORK array, and no error
64 ! message related to LWORK is issued by XERBLA.
65 !
66 ! INFO (output) INTEGER
67 ! = 0: successful exit
68 ! < 0: if INFO = -i, the i-th argument had an illegal value
69 ! > 0: if INFO = i, the algorithm failed to converge; i
70 ! off-diagonal elements of an intermediate tridiagonal
71 ! form did not converge to zero.
72 !
73 ! =====================================================================
74 !
75 ! .. Parameters ..
76 DOUBLE PRECISION ZERO, ONE
77 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
78 ! ..
79 ! .. Local Scalars ..
80 LOGICAL LOWER, LQUERY, WANTZ
81 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, &
82 LLWORK, LWKOPT, NB
83 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, &
84 SMLNUM
85 ! ..
86 ! .. External Functions ..
87 ! LOGICAL LSAME
88 ! INTEGER ILAENV
89 ! DOUBLE PRECISION DLAMCH, DLANSY
90 ! EXTERNAL LSAME
91 ! EXTERNAL 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 SUBROUTINE DSYEV