Annotation of rpl/lapack/lapack/dsyevd.f, revision 1.1.1.1
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsyevd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack