1: *> \brief <b> DSYEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSYEV + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyev.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyev.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyev.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER JOBZ, UPLO
25: * INTEGER INFO, LDA, LWORK, N
26: * ..
27: * .. Array Arguments ..
28: * DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
29: * ..
30: *
31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> DSYEV computes all eigenvalues and, optionally, eigenvectors of a
38: *> real symmetric matrix A.
39: *> \endverbatim
40: *
41: * Arguments:
42: * ==========
43: *
44: *> \param[in] JOBZ
45: *> \verbatim
46: *> JOBZ is CHARACTER*1
47: *> = 'N': Compute eigenvalues only;
48: *> = 'V': Compute eigenvalues and eigenvectors.
49: *> \endverbatim
50: *>
51: *> \param[in] UPLO
52: *> \verbatim
53: *> UPLO is CHARACTER*1
54: *> = 'U': Upper triangle of A is stored;
55: *> = 'L': Lower triangle of A is stored.
56: *> \endverbatim
57: *>
58: *> \param[in] N
59: *> \verbatim
60: *> N is INTEGER
61: *> The order of the matrix A. N >= 0.
62: *> \endverbatim
63: *>
64: *> \param[in,out] A
65: *> \verbatim
66: *> A is DOUBLE PRECISION array, dimension (LDA, N)
67: *> On entry, the symmetric matrix A. If UPLO = 'U', the
68: *> leading N-by-N upper triangular part of A contains the
69: *> upper triangular part of the matrix A. If UPLO = 'L',
70: *> the leading N-by-N lower triangular part of A contains
71: *> the lower triangular part of the matrix A.
72: *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
73: *> orthonormal eigenvectors of the matrix A.
74: *> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
75: *> or the upper triangle (if UPLO='U') of A, including the
76: *> diagonal, is destroyed.
77: *> \endverbatim
78: *>
79: *> \param[in] LDA
80: *> \verbatim
81: *> LDA is INTEGER
82: *> The leading dimension of the array A. LDA >= max(1,N).
83: *> \endverbatim
84: *>
85: *> \param[out] W
86: *> \verbatim
87: *> W is DOUBLE PRECISION array, dimension (N)
88: *> If INFO = 0, the eigenvalues in ascending order.
89: *> \endverbatim
90: *>
91: *> \param[out] WORK
92: *> \verbatim
93: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
94: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
95: *> \endverbatim
96: *>
97: *> \param[in] LWORK
98: *> \verbatim
99: *> LWORK is INTEGER
100: *> The length of the array WORK. LWORK >= max(1,3*N-1).
101: *> For optimal efficiency, LWORK >= (NB+2)*N,
102: *> where NB is the blocksize for DSYTRD returned by ILAENV.
103: *>
104: *> If LWORK = -1, then a workspace query is assumed; the routine
105: *> only calculates the optimal size of the WORK array, returns
106: *> this value as the first entry of the WORK array, and no error
107: *> message related to LWORK is issued by XERBLA.
108: *> \endverbatim
109: *>
110: *> \param[out] INFO
111: *> \verbatim
112: *> INFO is INTEGER
113: *> = 0: successful exit
114: *> < 0: if INFO = -i, the i-th argument had an illegal value
115: *> > 0: if INFO = i, the algorithm failed to converge; i
116: *> off-diagonal elements of an intermediate tridiagonal
117: *> form did not converge to zero.
118: *> \endverbatim
119: *
120: * Authors:
121: * ========
122: *
123: *> \author Univ. of Tennessee
124: *> \author Univ. of California Berkeley
125: *> \author Univ. of Colorado Denver
126: *> \author NAG Ltd.
127: *
128: *> \ingroup doubleSYeigen
129: *
130: * =====================================================================
131: SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
132: *
133: * -- LAPACK driver routine --
134: * -- LAPACK is a software package provided by Univ. of Tennessee, --
135: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136: *
137: * .. Scalar Arguments ..
138: CHARACTER JOBZ, UPLO
139: INTEGER INFO, LDA, LWORK, N
140: * ..
141: * .. Array Arguments ..
142: DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
143: * ..
144: *
145: * =====================================================================
146: *
147: * .. Parameters ..
148: DOUBLE PRECISION ZERO, ONE
149: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
150: * ..
151: * .. Local Scalars ..
152: LOGICAL LOWER, LQUERY, WANTZ
153: INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
154: $ LLWORK, LWKOPT, NB
155: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
156: $ SMLNUM
157: * ..
158: * .. External Functions ..
159: LOGICAL LSAME
160: INTEGER ILAENV
161: DOUBLE PRECISION DLAMCH, DLANSY
162: EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY
163: * ..
164: * .. External Subroutines ..
165: EXTERNAL DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD,
166: $ XERBLA
167: * ..
168: * .. Intrinsic Functions ..
169: INTRINSIC MAX, SQRT
170: * ..
171: * .. Executable Statements ..
172: *
173: * Test the input parameters.
174: *
175: WANTZ = LSAME( JOBZ, 'V' )
176: LOWER = LSAME( UPLO, 'L' )
177: LQUERY = ( LWORK.EQ.-1 )
178: *
179: INFO = 0
180: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
181: INFO = -1
182: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
183: INFO = -2
184: ELSE IF( N.LT.0 ) THEN
185: INFO = -3
186: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
187: INFO = -5
188: END IF
189: *
190: IF( INFO.EQ.0 ) THEN
191: NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
192: LWKOPT = MAX( 1, ( NB+2 )*N )
193: WORK( 1 ) = LWKOPT
194: *
195: IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY )
196: $ INFO = -8
197: END IF
198: *
199: IF( INFO.NE.0 ) THEN
200: CALL XERBLA( 'DSYEV ', -INFO )
201: RETURN
202: ELSE IF( LQUERY ) THEN
203: RETURN
204: END IF
205: *
206: * Quick return if possible
207: *
208: IF( N.EQ.0 ) THEN
209: RETURN
210: END IF
211: *
212: IF( N.EQ.1 ) THEN
213: W( 1 ) = A( 1, 1 )
214: WORK( 1 ) = 2
215: IF( WANTZ )
216: $ A( 1, 1 ) = ONE
217: RETURN
218: END IF
219: *
220: * Get machine constants.
221: *
222: SAFMIN = DLAMCH( 'Safe minimum' )
223: EPS = DLAMCH( 'Precision' )
224: SMLNUM = SAFMIN / EPS
225: BIGNUM = ONE / SMLNUM
226: RMIN = SQRT( SMLNUM )
227: RMAX = SQRT( BIGNUM )
228: *
229: * Scale matrix to allowable range, if necessary.
230: *
231: ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
232: ISCALE = 0
233: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
234: ISCALE = 1
235: SIGMA = RMIN / ANRM
236: ELSE IF( ANRM.GT.RMAX ) THEN
237: ISCALE = 1
238: SIGMA = RMAX / ANRM
239: END IF
240: IF( ISCALE.EQ.1 )
241: $ CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
242: *
243: * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
244: *
245: INDE = 1
246: INDTAU = INDE + N
247: INDWRK = INDTAU + N
248: LLWORK = LWORK - INDWRK + 1
249: CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
250: $ WORK( INDWRK ), LLWORK, IINFO )
251: *
252: * For eigenvalues only, call DSTERF. For eigenvectors, first call
253: * DORGTR to generate the orthogonal matrix, then call DSTEQR.
254: *
255: IF( .NOT.WANTZ ) THEN
256: CALL DSTERF( N, W, WORK( INDE ), INFO )
257: ELSE
258: CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
259: $ LLWORK, IINFO )
260: CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ),
261: $ INFO )
262: END IF
263: *
264: * If matrix was scaled, then rescale eigenvalues appropriately.
265: *
266: IF( ISCALE.EQ.1 ) THEN
267: IF( INFO.EQ.0 ) THEN
268: IMAX = N
269: ELSE
270: IMAX = INFO - 1
271: END IF
272: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
273: END IF
274: *
275: * Set WORK(1) to optimal workspace size.
276: *
277: WORK( 1 ) = LWKOPT
278: *
279: RETURN
280: *
281: * End of DSYEV
282: *
283: END
CVSweb interface <joel.bertrand@systella.fr>