Annotation of rpl/lapack/lapack/dstevx.f, revision 1.17
1.8 bertrand 1: *> \brief <b> DSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.15 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: *> \htmlonly
1.15 bertrand 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">
1.8 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.8 bertrand 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 )
1.15 bertrand 23: *
1.8 bertrand 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: * ..
1.15 bertrand 33: *
1.8 bertrand 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
1.13 bertrand 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'.
1.8 bertrand 95: *> \endverbatim
96: *>
97: *> \param[in] VU
98: *> \verbatim
99: *> VU is DOUBLE PRECISION
1.13 bertrand 100: *> If RANGE='V', the upper bound of the interval to
1.8 bertrand 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
1.13 bertrand 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'.
1.8 bertrand 112: *> \endverbatim
113: *>
114: *> \param[in] IU
115: *> \verbatim
116: *> IU is INTEGER
1.13 bertrand 117: *> If RANGE='I', the index of the
118: *> largest eigenvalue to be returned.
1.8 bertrand 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: *
1.15 bertrand 217: *> \author Univ. of Tennessee
218: *> \author Univ. of California Berkeley
219: *> \author Univ. of Colorado Denver
220: *> \author NAG Ltd.
1.8 bertrand 221: *
1.13 bertrand 222: *> \date June 2016
1.8 bertrand 223: *
224: *> \ingroup doubleOTHEReigen
225: *
226: * =====================================================================
1.1 bertrand 227: SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
228: $ M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO )
229: *
1.15 bertrand 230: * -- LAPACK driver routine (version 3.7.0) --
1.1 bertrand 231: * -- LAPACK is a software package provided by Univ. of Tennessee, --
232: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.13 bertrand 233: * June 2016
1.1 bertrand 234: *
235: * .. Scalar Arguments ..
236: CHARACTER JOBZ, RANGE
237: INTEGER IL, INFO, IU, LDZ, M, N
238: DOUBLE PRECISION ABSTOL, VL, VU
239: * ..
240: * .. Array Arguments ..
241: INTEGER IFAIL( * ), IWORK( * )
242: DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
243: * ..
244: *
245: * =====================================================================
246: *
247: * .. Parameters ..
248: DOUBLE PRECISION ZERO, ONE
249: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
250: * ..
251: * .. Local Scalars ..
252: LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
253: CHARACTER ORDER
254: INTEGER I, IMAX, INDIBL, INDISP, INDIWO, INDWRK,
255: $ ISCALE, ITMP1, J, JJ, NSPLIT
256: DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
257: $ TMP1, TNRM, VLL, VUU
258: * ..
259: * .. External Functions ..
260: LOGICAL LSAME
261: DOUBLE PRECISION DLAMCH, DLANST
262: EXTERNAL LSAME, DLAMCH, DLANST
263: * ..
264: * .. External Subroutines ..
265: EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTEIN, DSTEQR, DSTERF,
266: $ DSWAP, XERBLA
267: * ..
268: * .. Intrinsic Functions ..
269: INTRINSIC MAX, MIN, SQRT
270: * ..
271: * .. Executable Statements ..
272: *
273: * Test the input parameters.
274: *
275: WANTZ = LSAME( JOBZ, 'V' )
276: ALLEIG = LSAME( RANGE, 'A' )
277: VALEIG = LSAME( RANGE, 'V' )
278: INDEIG = LSAME( RANGE, 'I' )
279: *
280: INFO = 0
281: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
282: INFO = -1
283: ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
284: INFO = -2
285: ELSE IF( N.LT.0 ) THEN
286: INFO = -3
287: ELSE
288: IF( VALEIG ) THEN
289: IF( N.GT.0 .AND. VU.LE.VL )
290: $ INFO = -7
291: ELSE IF( INDEIG ) THEN
292: IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
293: INFO = -8
294: ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
295: INFO = -9
296: END IF
297: END IF
298: END IF
299: IF( INFO.EQ.0 ) THEN
300: IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) )
301: $ INFO = -14
302: END IF
303: *
304: IF( INFO.NE.0 ) THEN
305: CALL XERBLA( 'DSTEVX', -INFO )
306: RETURN
307: END IF
308: *
309: * Quick return if possible
310: *
311: M = 0
312: IF( N.EQ.0 )
313: $ RETURN
314: *
315: IF( N.EQ.1 ) THEN
316: IF( ALLEIG .OR. INDEIG ) THEN
317: M = 1
318: W( 1 ) = D( 1 )
319: ELSE
320: IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN
321: M = 1
322: W( 1 ) = D( 1 )
323: END IF
324: END IF
325: IF( WANTZ )
326: $ Z( 1, 1 ) = ONE
327: RETURN
328: END IF
329: *
330: * Get machine constants.
331: *
332: SAFMIN = DLAMCH( 'Safe minimum' )
333: EPS = DLAMCH( 'Precision' )
334: SMLNUM = SAFMIN / EPS
335: BIGNUM = ONE / SMLNUM
336: RMIN = SQRT( SMLNUM )
337: RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
338: *
339: * Scale matrix to allowable range, if necessary.
340: *
341: ISCALE = 0
342: IF( VALEIG ) THEN
343: VLL = VL
344: VUU = VU
345: ELSE
346: VLL = ZERO
347: VUU = ZERO
348: END IF
349: TNRM = DLANST( 'M', N, D, E )
350: IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
351: ISCALE = 1
352: SIGMA = RMIN / TNRM
353: ELSE IF( TNRM.GT.RMAX ) THEN
354: ISCALE = 1
355: SIGMA = RMAX / TNRM
356: END IF
357: IF( ISCALE.EQ.1 ) THEN
358: CALL DSCAL( N, SIGMA, D, 1 )
359: CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
360: IF( VALEIG ) THEN
361: VLL = VL*SIGMA
362: VUU = VU*SIGMA
363: END IF
364: END IF
365: *
366: * If all eigenvalues are desired and ABSTOL is less than zero, then
367: * call DSTERF or SSTEQR. If this fails for some eigenvalue, then
368: * try DSTEBZ.
369: *
370: TEST = .FALSE.
371: IF( INDEIG ) THEN
372: IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
373: TEST = .TRUE.
374: END IF
375: END IF
376: IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN
377: CALL DCOPY( N, D, 1, W, 1 )
378: CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 )
379: INDWRK = N + 1
380: IF( .NOT.WANTZ ) THEN
381: CALL DSTERF( N, W, WORK, INFO )
382: ELSE
383: CALL DSTEQR( 'I', N, W, WORK, Z, LDZ, WORK( INDWRK ), INFO )
384: IF( INFO.EQ.0 ) THEN
385: DO 10 I = 1, N
386: IFAIL( I ) = 0
387: 10 CONTINUE
388: END IF
389: END IF
390: IF( INFO.EQ.0 ) THEN
391: M = N
392: GO TO 20
393: END IF
394: INFO = 0
395: END IF
396: *
397: * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
398: *
399: IF( WANTZ ) THEN
400: ORDER = 'B'
401: ELSE
402: ORDER = 'E'
403: END IF
404: INDWRK = 1
405: INDIBL = 1
406: INDISP = INDIBL + N
407: INDIWO = INDISP + N
408: CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M,
409: $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ),
410: $ WORK( INDWRK ), IWORK( INDIWO ), INFO )
411: *
412: IF( WANTZ ) THEN
413: CALL DSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ),
414: $ Z, LDZ, WORK( INDWRK ), IWORK( INDIWO ), IFAIL,
415: $ INFO )
416: END IF
417: *
418: * If matrix was scaled, then rescale eigenvalues appropriately.
419: *
420: 20 CONTINUE
421: IF( ISCALE.EQ.1 ) THEN
422: IF( INFO.EQ.0 ) THEN
423: IMAX = M
424: ELSE
425: IMAX = INFO - 1
426: END IF
427: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
428: END IF
429: *
430: * If eigenvalues are not in order, then sort them, along with
431: * eigenvectors.
432: *
433: IF( WANTZ ) THEN
434: DO 40 J = 1, M - 1
435: I = 0
436: TMP1 = W( J )
437: DO 30 JJ = J + 1, M
438: IF( W( JJ ).LT.TMP1 ) THEN
439: I = JJ
440: TMP1 = W( JJ )
441: END IF
442: 30 CONTINUE
443: *
444: IF( I.NE.0 ) THEN
445: ITMP1 = IWORK( INDIBL+I-1 )
446: W( I ) = W( J )
447: IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
448: W( J ) = TMP1
449: IWORK( INDIBL+J-1 ) = ITMP1
450: CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
451: IF( INFO.NE.0 ) THEN
452: ITMP1 = IFAIL( I )
453: IFAIL( I ) = IFAIL( J )
454: IFAIL( J ) = ITMP1
455: END IF
456: END IF
457: 40 CONTINUE
458: END IF
459: *
460: RETURN
461: *
462: * End of DSTEVX
463: *
464: 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