1: *> \brief <b> DSTEVX 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 DSTEVX + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstevx.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstevx.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstevx.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
22: * M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER JOBZ, RANGE
26: * INTEGER IL, INFO, IU, LDZ, M, N
27: * DOUBLE PRECISION ABSTOL, VL, VU
28: * ..
29: * .. Array Arguments ..
30: * INTEGER IFAIL( * ), IWORK( * )
31: * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DSTEVX computes selected eigenvalues and, optionally, eigenvectors
41: *> of a real symmetric tridiagonal matrix A. Eigenvalues and
42: *> eigenvectors can be selected by specifying either a range of values
43: *> or a range of indices for the desired eigenvalues.
44: *> \endverbatim
45: *
46: * Arguments:
47: * ==========
48: *
49: *> \param[in] JOBZ
50: *> \verbatim
51: *> JOBZ is CHARACTER*1
52: *> = 'N': Compute eigenvalues only;
53: *> = 'V': Compute eigenvalues and eigenvectors.
54: *> \endverbatim
55: *>
56: *> \param[in] RANGE
57: *> \verbatim
58: *> RANGE is CHARACTER*1
59: *> = 'A': all eigenvalues will be found.
60: *> = 'V': all eigenvalues in the half-open interval (VL,VU]
61: *> will be found.
62: *> = 'I': the IL-th through IU-th eigenvalues will be found.
63: *> \endverbatim
64: *>
65: *> \param[in] N
66: *> \verbatim
67: *> N is INTEGER
68: *> The order of the matrix. N >= 0.
69: *> \endverbatim
70: *>
71: *> \param[in,out] D
72: *> \verbatim
73: *> D is DOUBLE PRECISION array, dimension (N)
74: *> On entry, the n diagonal elements of the tridiagonal matrix
75: *> A.
76: *> On exit, D may be multiplied by a constant factor chosen
77: *> to avoid over/underflow in computing the eigenvalues.
78: *> \endverbatim
79: *>
80: *> \param[in,out] E
81: *> \verbatim
82: *> E is DOUBLE PRECISION array, dimension (max(1,N-1))
83: *> On entry, the (n-1) subdiagonal elements of the tridiagonal
84: *> matrix A in elements 1 to N-1 of E.
85: *> On exit, E may be multiplied by a constant factor chosen
86: *> to avoid over/underflow in computing the eigenvalues.
87: *> \endverbatim
88: *>
89: *> \param[in] VL
90: *> \verbatim
91: *> VL is DOUBLE PRECISION
92: *> If RANGE='V', the lower bound of the interval to
93: *> be searched for eigenvalues. VL < VU.
94: *> Not referenced if RANGE = 'A' or 'I'.
95: *> \endverbatim
96: *>
97: *> \param[in] VU
98: *> \verbatim
99: *> VU is DOUBLE PRECISION
100: *> If RANGE='V', the upper bound of the interval to
101: *> be searched for eigenvalues. VL < VU.
102: *> Not referenced if RANGE = 'A' or 'I'.
103: *> \endverbatim
104: *>
105: *> \param[in] IL
106: *> \verbatim
107: *> IL is INTEGER
108: *> If RANGE='I', the index of the
109: *> smallest eigenvalue to be returned.
110: *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
111: *> Not referenced if RANGE = 'A' or 'V'.
112: *> \endverbatim
113: *>
114: *> \param[in] IU
115: *> \verbatim
116: *> IU is INTEGER
117: *> If RANGE='I', the index of the
118: *> largest eigenvalue to be returned.
119: *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
120: *> Not referenced if RANGE = 'A' or 'V'.
121: *> \endverbatim
122: *>
123: *> \param[in] ABSTOL
124: *> \verbatim
125: *> ABSTOL is DOUBLE PRECISION
126: *> The absolute error tolerance for the eigenvalues.
127: *> An approximate eigenvalue is accepted as converged
128: *> when it is determined to lie in an interval [a,b]
129: *> of width less than or equal to
130: *>
131: *> ABSTOL + EPS * max( |a|,|b| ) ,
132: *>
133: *> where EPS is the machine precision. If ABSTOL is less
134: *> than or equal to zero, then EPS*|T| will be used in
135: *> its place, where |T| is the 1-norm of the tridiagonal
136: *> matrix.
137: *>
138: *> Eigenvalues will be computed most accurately when ABSTOL is
139: *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
140: *> If this routine returns with INFO>0, indicating that some
141: *> eigenvectors did not converge, try setting ABSTOL to
142: *> 2*DLAMCH('S').
143: *>
144: *> See "Computing Small Singular Values of Bidiagonal Matrices
145: *> with Guaranteed High Relative Accuracy," by Demmel and
146: *> Kahan, LAPACK Working Note #3.
147: *> \endverbatim
148: *>
149: *> \param[out] M
150: *> \verbatim
151: *> M is INTEGER
152: *> The total number of eigenvalues found. 0 <= M <= N.
153: *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
154: *> \endverbatim
155: *>
156: *> \param[out] W
157: *> \verbatim
158: *> W is DOUBLE PRECISION array, dimension (N)
159: *> The first M elements contain the selected eigenvalues in
160: *> ascending order.
161: *> \endverbatim
162: *>
163: *> \param[out] Z
164: *> \verbatim
165: *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
166: *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
167: *> contain the orthonormal eigenvectors of the matrix A
168: *> corresponding to the selected eigenvalues, with the i-th
169: *> column of Z holding the eigenvector associated with W(i).
170: *> If an eigenvector fails to converge (INFO > 0), then that
171: *> column of Z contains the latest approximation to the
172: *> eigenvector, and the index of the eigenvector is returned
173: *> in IFAIL. If JOBZ = 'N', then Z is not referenced.
174: *> Note: the user must ensure that at least max(1,M) columns are
175: *> supplied in the array Z; if RANGE = 'V', the exact value of M
176: *> is not known in advance and an upper bound must be used.
177: *> \endverbatim
178: *>
179: *> \param[in] LDZ
180: *> \verbatim
181: *> LDZ is INTEGER
182: *> The leading dimension of the array Z. LDZ >= 1, and if
183: *> JOBZ = 'V', LDZ >= max(1,N).
184: *> \endverbatim
185: *>
186: *> \param[out] WORK
187: *> \verbatim
188: *> WORK is DOUBLE PRECISION array, dimension (5*N)
189: *> \endverbatim
190: *>
191: *> \param[out] IWORK
192: *> \verbatim
193: *> IWORK is INTEGER array, dimension (5*N)
194: *> \endverbatim
195: *>
196: *> \param[out] IFAIL
197: *> \verbatim
198: *> IFAIL is INTEGER array, dimension (N)
199: *> If JOBZ = 'V', then if INFO = 0, the first M elements of
200: *> IFAIL are zero. If INFO > 0, then IFAIL contains the
201: *> indices of the eigenvectors that failed to converge.
202: *> If JOBZ = 'N', then IFAIL is not referenced.
203: *> \endverbatim
204: *>
205: *> \param[out] INFO
206: *> \verbatim
207: *> INFO is INTEGER
208: *> = 0: successful exit
209: *> < 0: if INFO = -i, the i-th argument had an illegal value
210: *> > 0: if INFO = i, then i eigenvectors failed to converge.
211: *> Their indices are stored in array IFAIL.
212: *> \endverbatim
213: *
214: * Authors:
215: * ========
216: *
217: *> \author Univ. of Tennessee
218: *> \author Univ. of California Berkeley
219: *> \author Univ. of Colorado Denver
220: *> \author NAG Ltd.
221: *
222: *> \ingroup doubleOTHEReigen
223: *
224: * =====================================================================
225: SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
226: $ M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO )
227: *
228: * -- LAPACK driver routine --
229: * -- LAPACK is a software package provided by Univ. of Tennessee, --
230: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231: *
232: * .. Scalar Arguments ..
233: CHARACTER JOBZ, RANGE
234: INTEGER IL, INFO, IU, LDZ, M, N
235: DOUBLE PRECISION ABSTOL, VL, VU
236: * ..
237: * .. Array Arguments ..
238: INTEGER IFAIL( * ), IWORK( * )
239: DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
240: * ..
241: *
242: * =====================================================================
243: *
244: * .. Parameters ..
245: DOUBLE PRECISION ZERO, ONE
246: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
247: * ..
248: * .. Local Scalars ..
249: LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
250: CHARACTER ORDER
251: INTEGER I, IMAX, INDIBL, INDISP, INDIWO, INDWRK,
252: $ ISCALE, ITMP1, J, JJ, NSPLIT
253: DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
254: $ TMP1, TNRM, VLL, VUU
255: * ..
256: * .. External Functions ..
257: LOGICAL LSAME
258: DOUBLE PRECISION DLAMCH, DLANST
259: EXTERNAL LSAME, DLAMCH, DLANST
260: * ..
261: * .. External Subroutines ..
262: EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTEIN, DSTEQR, DSTERF,
263: $ DSWAP, XERBLA
264: * ..
265: * .. Intrinsic Functions ..
266: INTRINSIC MAX, MIN, SQRT
267: * ..
268: * .. Executable Statements ..
269: *
270: * Test the input parameters.
271: *
272: WANTZ = LSAME( JOBZ, 'V' )
273: ALLEIG = LSAME( RANGE, 'A' )
274: VALEIG = LSAME( RANGE, 'V' )
275: INDEIG = LSAME( RANGE, 'I' )
276: *
277: INFO = 0
278: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
279: INFO = -1
280: ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
281: INFO = -2
282: ELSE IF( N.LT.0 ) THEN
283: INFO = -3
284: ELSE
285: IF( VALEIG ) THEN
286: IF( N.GT.0 .AND. VU.LE.VL )
287: $ INFO = -7
288: ELSE IF( INDEIG ) THEN
289: IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
290: INFO = -8
291: ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
292: INFO = -9
293: END IF
294: END IF
295: END IF
296: IF( INFO.EQ.0 ) THEN
297: IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) )
298: $ INFO = -14
299: END IF
300: *
301: IF( INFO.NE.0 ) THEN
302: CALL XERBLA( 'DSTEVX', -INFO )
303: RETURN
304: END IF
305: *
306: * Quick return if possible
307: *
308: M = 0
309: IF( N.EQ.0 )
310: $ RETURN
311: *
312: IF( N.EQ.1 ) THEN
313: IF( ALLEIG .OR. INDEIG ) THEN
314: M = 1
315: W( 1 ) = D( 1 )
316: ELSE
317: IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN
318: M = 1
319: W( 1 ) = D( 1 )
320: END IF
321: END IF
322: IF( WANTZ )
323: $ Z( 1, 1 ) = ONE
324: RETURN
325: END IF
326: *
327: * Get machine constants.
328: *
329: SAFMIN = DLAMCH( 'Safe minimum' )
330: EPS = DLAMCH( 'Precision' )
331: SMLNUM = SAFMIN / EPS
332: BIGNUM = ONE / SMLNUM
333: RMIN = SQRT( SMLNUM )
334: RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
335: *
336: * Scale matrix to allowable range, if necessary.
337: *
338: ISCALE = 0
339: IF( VALEIG ) THEN
340: VLL = VL
341: VUU = VU
342: ELSE
343: VLL = ZERO
344: VUU = ZERO
345: END IF
346: TNRM = DLANST( 'M', N, D, E )
347: IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
348: ISCALE = 1
349: SIGMA = RMIN / TNRM
350: ELSE IF( TNRM.GT.RMAX ) THEN
351: ISCALE = 1
352: SIGMA = RMAX / TNRM
353: END IF
354: IF( ISCALE.EQ.1 ) THEN
355: CALL DSCAL( N, SIGMA, D, 1 )
356: CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
357: IF( VALEIG ) THEN
358: VLL = VL*SIGMA
359: VUU = VU*SIGMA
360: END IF
361: END IF
362: *
363: * If all eigenvalues are desired and ABSTOL is less than zero, then
364: * call DSTERF or SSTEQR. If this fails for some eigenvalue, then
365: * try DSTEBZ.
366: *
367: TEST = .FALSE.
368: IF( INDEIG ) THEN
369: IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
370: TEST = .TRUE.
371: END IF
372: END IF
373: IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN
374: CALL DCOPY( N, D, 1, W, 1 )
375: CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 )
376: INDWRK = N + 1
377: IF( .NOT.WANTZ ) THEN
378: CALL DSTERF( N, W, WORK, INFO )
379: ELSE
380: CALL DSTEQR( 'I', N, W, WORK, Z, LDZ, WORK( INDWRK ), INFO )
381: IF( INFO.EQ.0 ) THEN
382: DO 10 I = 1, N
383: IFAIL( I ) = 0
384: 10 CONTINUE
385: END IF
386: END IF
387: IF( INFO.EQ.0 ) THEN
388: M = N
389: GO TO 20
390: END IF
391: INFO = 0
392: END IF
393: *
394: * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
395: *
396: IF( WANTZ ) THEN
397: ORDER = 'B'
398: ELSE
399: ORDER = 'E'
400: END IF
401: INDWRK = 1
402: INDIBL = 1
403: INDISP = INDIBL + N
404: INDIWO = INDISP + N
405: CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M,
406: $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ),
407: $ WORK( INDWRK ), IWORK( INDIWO ), INFO )
408: *
409: IF( WANTZ ) THEN
410: CALL DSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ),
411: $ Z, LDZ, WORK( INDWRK ), IWORK( INDIWO ), IFAIL,
412: $ INFO )
413: END IF
414: *
415: * If matrix was scaled, then rescale eigenvalues appropriately.
416: *
417: 20 CONTINUE
418: IF( ISCALE.EQ.1 ) THEN
419: IF( INFO.EQ.0 ) THEN
420: IMAX = M
421: ELSE
422: IMAX = INFO - 1
423: END IF
424: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
425: END IF
426: *
427: * If eigenvalues are not in order, then sort them, along with
428: * eigenvectors.
429: *
430: IF( WANTZ ) THEN
431: DO 40 J = 1, M - 1
432: I = 0
433: TMP1 = W( J )
434: DO 30 JJ = J + 1, M
435: IF( W( JJ ).LT.TMP1 ) THEN
436: I = JJ
437: TMP1 = W( JJ )
438: END IF
439: 30 CONTINUE
440: *
441: IF( I.NE.0 ) THEN
442: ITMP1 = IWORK( INDIBL+I-1 )
443: W( I ) = W( J )
444: IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
445: W( J ) = TMP1
446: IWORK( INDIBL+J-1 ) = ITMP1
447: CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
448: IF( INFO.NE.0 ) THEN
449: ITMP1 = IFAIL( I )
450: IFAIL( I ) = IFAIL( J )
451: IFAIL( J ) = ITMP1
452: END IF
453: END IF
454: 40 CONTINUE
455: END IF
456: *
457: RETURN
458: *
459: * End of DSTEVX
460: *
461: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dstevx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack