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: *> \date November 2011
129: *
130: *> \ingroup doubleSYeigen
131: *
132: * =====================================================================
133: SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
134: *
135: * -- LAPACK driver routine (version 3.4.0) --
136: * -- LAPACK is a software package provided by Univ. of Tennessee, --
137: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138: * November 2011
139: *
140: * .. Scalar Arguments ..
141: CHARACTER JOBZ, UPLO
142: INTEGER INFO, LDA, LWORK, N
143: * ..
144: * .. Array Arguments ..
145: DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
146: * ..
147: *
148: * =====================================================================
149: *
150: * .. Parameters ..
151: DOUBLE PRECISION ZERO, ONE
152: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
153: * ..
154: * .. Local Scalars ..
155: LOGICAL LOWER, LQUERY, WANTZ
156: INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
157: $ LLWORK, LWKOPT, NB
158: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
159: $ SMLNUM
160: * ..
161: * .. External Functions ..
162: LOGICAL LSAME
163: INTEGER ILAENV
164: DOUBLE PRECISION DLAMCH, DLANSY
165: EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY
166: * ..
167: * .. External Subroutines ..
168: EXTERNAL DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD,
169: $ XERBLA
170: * ..
171: * .. Intrinsic Functions ..
172: INTRINSIC MAX, SQRT
173: * ..
174: * .. Executable Statements ..
175: *
176: * Test the input parameters.
177: *
178: WANTZ = LSAME( JOBZ, 'V' )
179: LOWER = LSAME( UPLO, 'L' )
180: LQUERY = ( LWORK.EQ.-1 )
181: *
182: INFO = 0
183: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
184: INFO = -1
185: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
186: INFO = -2
187: ELSE IF( N.LT.0 ) THEN
188: INFO = -3
189: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
190: INFO = -5
191: END IF
192: *
193: IF( INFO.EQ.0 ) THEN
194: NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
195: LWKOPT = MAX( 1, ( NB+2 )*N )
196: WORK( 1 ) = LWKOPT
197: *
198: IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY )
199: $ INFO = -8
200: END IF
201: *
202: IF( INFO.NE.0 ) THEN
203: CALL XERBLA( 'DSYEV ', -INFO )
204: RETURN
205: ELSE IF( LQUERY ) THEN
206: RETURN
207: END IF
208: *
209: * Quick return if possible
210: *
211: IF( N.EQ.0 ) THEN
212: RETURN
213: END IF
214: *
215: IF( N.EQ.1 ) THEN
216: W( 1 ) = A( 1, 1 )
217: WORK( 1 ) = 2
218: IF( WANTZ )
219: $ A( 1, 1 ) = ONE
220: RETURN
221: END IF
222: *
223: * Get machine constants.
224: *
225: SAFMIN = DLAMCH( 'Safe minimum' )
226: EPS = DLAMCH( 'Precision' )
227: SMLNUM = SAFMIN / EPS
228: BIGNUM = ONE / SMLNUM
229: RMIN = SQRT( SMLNUM )
230: RMAX = SQRT( BIGNUM )
231: *
232: * Scale matrix to allowable range, if necessary.
233: *
234: ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
235: ISCALE = 0
236: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
237: ISCALE = 1
238: SIGMA = RMIN / ANRM
239: ELSE IF( ANRM.GT.RMAX ) THEN
240: ISCALE = 1
241: SIGMA = RMAX / ANRM
242: END IF
243: IF( ISCALE.EQ.1 )
244: $ CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
245: *
246: * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
247: *
248: INDE = 1
249: INDTAU = INDE + N
250: INDWRK = INDTAU + N
251: LLWORK = LWORK - INDWRK + 1
252: CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
253: $ WORK( INDWRK ), LLWORK, IINFO )
254: *
255: * For eigenvalues only, call DSTERF. For eigenvectors, first call
256: * DORGTR to generate the orthogonal matrix, then call DSTEQR.
257: *
258: IF( .NOT.WANTZ ) THEN
259: CALL DSTERF( N, W, WORK( INDE ), INFO )
260: ELSE
261: CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
262: $ LLWORK, IINFO )
263: CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ),
264: $ INFO )
265: END IF
266: *
267: * If matrix was scaled, then rescale eigenvalues appropriately.
268: *
269: IF( ISCALE.EQ.1 ) THEN
270: IF( INFO.EQ.0 ) THEN
271: IMAX = N
272: ELSE
273: IMAX = INFO - 1
274: END IF
275: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
276: END IF
277: *
278: * Set WORK(1) to optimal workspace size.
279: *
280: WORK( 1 ) = LWKOPT
281: *
282: RETURN
283: *
284: * End of DSYEV
285: *
286: END
CVSweb interface <joel.bertrand@systella.fr>