1: *> \brief <b> DSPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER 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 DSPEV + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dspev.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dspev.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dspev.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER JOBZ, UPLO
25: * INTEGER INFO, LDZ, N
26: * ..
27: * .. Array Arguments ..
28: * DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
29: * ..
30: *
31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> DSPEV computes all the eigenvalues and, optionally, eigenvectors of a
38: *> real symmetric matrix A in packed storage.
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] AP
65: *> \verbatim
66: *> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
67: *> On entry, the upper or lower triangle of the symmetric matrix
68: *> A, packed columnwise in a linear array. The j-th column of A
69: *> is stored in the array AP as follows:
70: *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
71: *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
72: *>
73: *> On exit, AP is overwritten by values generated during the
74: *> reduction to tridiagonal form. If UPLO = 'U', the diagonal
75: *> and first superdiagonal of the tridiagonal matrix T overwrite
76: *> the corresponding elements of A, and if UPLO = 'L', the
77: *> diagonal and first subdiagonal of T overwrite the
78: *> corresponding elements of A.
79: *> \endverbatim
80: *>
81: *> \param[out] W
82: *> \verbatim
83: *> W is DOUBLE PRECISION array, dimension (N)
84: *> If INFO = 0, the eigenvalues in ascending order.
85: *> \endverbatim
86: *>
87: *> \param[out] Z
88: *> \verbatim
89: *> Z is DOUBLE PRECISION array, dimension (LDZ, N)
90: *> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
91: *> eigenvectors of the matrix A, with the i-th column of Z
92: *> holding the eigenvector associated with W(i).
93: *> If JOBZ = 'N', then Z is not referenced.
94: *> \endverbatim
95: *>
96: *> \param[in] LDZ
97: *> \verbatim
98: *> LDZ is INTEGER
99: *> The leading dimension of the array Z. LDZ >= 1, and if
100: *> JOBZ = 'V', LDZ >= max(1,N).
101: *> \endverbatim
102: *>
103: *> \param[out] WORK
104: *> \verbatim
105: *> WORK is DOUBLE PRECISION array, dimension (3*N)
106: *> \endverbatim
107: *>
108: *> \param[out] INFO
109: *> \verbatim
110: *> INFO is INTEGER
111: *> = 0: successful exit.
112: *> < 0: if INFO = -i, the i-th argument had an illegal value.
113: *> > 0: if INFO = i, the algorithm failed to converge; i
114: *> off-diagonal elements of an intermediate tridiagonal
115: *> form did not converge to zero.
116: *> \endverbatim
117: *
118: * Authors:
119: * ========
120: *
121: *> \author Univ. of Tennessee
122: *> \author Univ. of California Berkeley
123: *> \author Univ. of Colorado Denver
124: *> \author NAG Ltd.
125: *
126: *> \ingroup doubleOTHEReigen
127: *
128: * =====================================================================
129: SUBROUTINE DSPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, INFO )
130: *
131: * -- LAPACK driver routine --
132: * -- LAPACK is a software package provided by Univ. of Tennessee, --
133: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134: *
135: * .. Scalar Arguments ..
136: CHARACTER JOBZ, UPLO
137: INTEGER INFO, LDZ, N
138: * ..
139: * .. Array Arguments ..
140: DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
141: * ..
142: *
143: * =====================================================================
144: *
145: * .. Parameters ..
146: DOUBLE PRECISION ZERO, ONE
147: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
148: * ..
149: * .. Local Scalars ..
150: LOGICAL WANTZ
151: INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE
152: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
153: $ SMLNUM
154: * ..
155: * .. External Functions ..
156: LOGICAL LSAME
157: DOUBLE PRECISION DLAMCH, DLANSP
158: EXTERNAL LSAME, DLAMCH, DLANSP
159: * ..
160: * .. External Subroutines ..
161: EXTERNAL DOPGTR, DSCAL, DSPTRD, DSTEQR, DSTERF, XERBLA
162: * ..
163: * .. Intrinsic Functions ..
164: INTRINSIC SQRT
165: * ..
166: * .. Executable Statements ..
167: *
168: * Test the input parameters.
169: *
170: WANTZ = LSAME( JOBZ, 'V' )
171: *
172: INFO = 0
173: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
174: INFO = -1
175: ELSE IF( .NOT.( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) )
176: $ THEN
177: INFO = -2
178: ELSE IF( N.LT.0 ) THEN
179: INFO = -3
180: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
181: INFO = -7
182: END IF
183: *
184: IF( INFO.NE.0 ) THEN
185: CALL XERBLA( 'DSPEV ', -INFO )
186: RETURN
187: END IF
188: *
189: * Quick return if possible
190: *
191: IF( N.EQ.0 )
192: $ RETURN
193: *
194: IF( N.EQ.1 ) THEN
195: W( 1 ) = AP( 1 )
196: IF( WANTZ )
197: $ Z( 1, 1 ) = ONE
198: RETURN
199: END IF
200: *
201: * Get machine constants.
202: *
203: SAFMIN = DLAMCH( 'Safe minimum' )
204: EPS = DLAMCH( 'Precision' )
205: SMLNUM = SAFMIN / EPS
206: BIGNUM = ONE / SMLNUM
207: RMIN = SQRT( SMLNUM )
208: RMAX = SQRT( BIGNUM )
209: *
210: * Scale matrix to allowable range, if necessary.
211: *
212: ANRM = DLANSP( 'M', UPLO, N, AP, WORK )
213: ISCALE = 0
214: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
215: ISCALE = 1
216: SIGMA = RMIN / ANRM
217: ELSE IF( ANRM.GT.RMAX ) THEN
218: ISCALE = 1
219: SIGMA = RMAX / ANRM
220: END IF
221: IF( ISCALE.EQ.1 ) THEN
222: CALL DSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
223: END IF
224: *
225: * Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
226: *
227: INDE = 1
228: INDTAU = INDE + N
229: CALL DSPTRD( UPLO, N, AP, W, WORK( INDE ), WORK( INDTAU ), IINFO )
230: *
231: * For eigenvalues only, call DSTERF. For eigenvectors, first call
232: * DOPGTR to generate the orthogonal matrix, then call DSTEQR.
233: *
234: IF( .NOT.WANTZ ) THEN
235: CALL DSTERF( N, W, WORK( INDE ), INFO )
236: ELSE
237: INDWRK = INDTAU + N
238: CALL DOPGTR( UPLO, N, AP, WORK( INDTAU ), Z, LDZ,
239: $ WORK( INDWRK ), IINFO )
240: CALL DSTEQR( JOBZ, N, W, WORK( INDE ), Z, LDZ, WORK( INDTAU ),
241: $ INFO )
242: END IF
243: *
244: * If matrix was scaled, then rescale eigenvalues appropriately.
245: *
246: IF( ISCALE.EQ.1 ) THEN
247: IF( INFO.EQ.0 ) THEN
248: IMAX = N
249: ELSE
250: IMAX = INFO - 1
251: END IF
252: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
253: END IF
254: *
255: RETURN
256: *
257: * End of DSPEV
258: *
259: END
CVSweb interface <joel.bertrand@systella.fr>