Annotation of rpl/lapack/lapack/dstevr.f, revision 1.8
1.8 ! bertrand 1: *> \brief <b> DSTEVR 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 DSTEVR + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstevr.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstevr.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstevr.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
! 22: * M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK,
! 23: * LIWORK, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER JOBZ, RANGE
! 27: * INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
! 28: * DOUBLE PRECISION ABSTOL, VL, VU
! 29: * ..
! 30: * .. Array Arguments ..
! 31: * INTEGER ISUPPZ( * ), IWORK( * )
! 32: * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
! 33: * ..
! 34: *
! 35: *
! 36: *> \par Purpose:
! 37: * =============
! 38: *>
! 39: *> \verbatim
! 40: *>
! 41: *> DSTEVR computes selected eigenvalues and, optionally, eigenvectors
! 42: *> of a real symmetric tridiagonal matrix T. Eigenvalues and
! 43: *> eigenvectors can be selected by specifying either a range of values
! 44: *> or a range of indices for the desired eigenvalues.
! 45: *>
! 46: *> Whenever possible, DSTEVR calls DSTEMR to compute the
! 47: *> eigenspectrum using Relatively Robust Representations. DSTEMR
! 48: *> computes eigenvalues by the dqds algorithm, while orthogonal
! 49: *> eigenvectors are computed from various "good" L D L^T representations
! 50: *> (also known as Relatively Robust Representations). Gram-Schmidt
! 51: *> orthogonalization is avoided as far as possible. More specifically,
! 52: *> the various steps of the algorithm are as follows. For the i-th
! 53: *> unreduced block of T,
! 54: *> (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
! 55: *> is a relatively robust representation,
! 56: *> (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
! 57: *> relative accuracy by the dqds algorithm,
! 58: *> (c) If there is a cluster of close eigenvalues, "choose" sigma_i
! 59: *> close to the cluster, and go to step (a),
! 60: *> (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
! 61: *> compute the corresponding eigenvector by forming a
! 62: *> rank-revealing twisted factorization.
! 63: *> The desired accuracy of the output can be specified by the input
! 64: *> parameter ABSTOL.
! 65: *>
! 66: *> For more details, see "A new O(n^2) algorithm for the symmetric
! 67: *> tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
! 68: *> Computer Science Division Technical Report No. UCB//CSD-97-971,
! 69: *> UC Berkeley, May 1997.
! 70: *>
! 71: *>
! 72: *> Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested
! 73: *> on machines which conform to the ieee-754 floating point standard.
! 74: *> DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
! 75: *> when partial spectrum requests are made.
! 76: *>
! 77: *> Normal execution of DSTEMR may create NaNs and infinities and
! 78: *> hence may abort due to a floating point exception in environments
! 79: *> which do not handle NaNs and infinities in the ieee standard default
! 80: *> manner.
! 81: *> \endverbatim
! 82: *
! 83: * Arguments:
! 84: * ==========
! 85: *
! 86: *> \param[in] JOBZ
! 87: *> \verbatim
! 88: *> JOBZ is CHARACTER*1
! 89: *> = 'N': Compute eigenvalues only;
! 90: *> = 'V': Compute eigenvalues and eigenvectors.
! 91: *> \endverbatim
! 92: *>
! 93: *> \param[in] RANGE
! 94: *> \verbatim
! 95: *> RANGE is CHARACTER*1
! 96: *> = 'A': all eigenvalues will be found.
! 97: *> = 'V': all eigenvalues in the half-open interval (VL,VU]
! 98: *> will be found.
! 99: *> = 'I': the IL-th through IU-th eigenvalues will be found.
! 100: *> For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
! 101: *> DSTEIN are called
! 102: *> \endverbatim
! 103: *>
! 104: *> \param[in] N
! 105: *> \verbatim
! 106: *> N is INTEGER
! 107: *> The order of the matrix. N >= 0.
! 108: *> \endverbatim
! 109: *>
! 110: *> \param[in,out] D
! 111: *> \verbatim
! 112: *> D is DOUBLE PRECISION array, dimension (N)
! 113: *> On entry, the n diagonal elements of the tridiagonal matrix
! 114: *> A.
! 115: *> On exit, D may be multiplied by a constant factor chosen
! 116: *> to avoid over/underflow in computing the eigenvalues.
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[in,out] E
! 120: *> \verbatim
! 121: *> E is DOUBLE PRECISION array, dimension (max(1,N-1))
! 122: *> On entry, the (n-1) subdiagonal elements of the tridiagonal
! 123: *> matrix A in elements 1 to N-1 of E.
! 124: *> On exit, E may be multiplied by a constant factor chosen
! 125: *> to avoid over/underflow in computing the eigenvalues.
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[in] VL
! 129: *> \verbatim
! 130: *> VL is DOUBLE PRECISION
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[in] VU
! 134: *> \verbatim
! 135: *> VU is DOUBLE PRECISION
! 136: *> If RANGE='V', the lower and upper bounds of the interval to
! 137: *> be searched for eigenvalues. VL < VU.
! 138: *> Not referenced if RANGE = 'A' or 'I'.
! 139: *> \endverbatim
! 140: *>
! 141: *> \param[in] IL
! 142: *> \verbatim
! 143: *> IL is INTEGER
! 144: *> \endverbatim
! 145: *>
! 146: *> \param[in] IU
! 147: *> \verbatim
! 148: *> IU is INTEGER
! 149: *> If RANGE='I', the indices (in ascending order) of the
! 150: *> smallest and largest eigenvalues to be returned.
! 151: *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
! 152: *> Not referenced if RANGE = 'A' or 'V'.
! 153: *> \endverbatim
! 154: *>
! 155: *> \param[in] ABSTOL
! 156: *> \verbatim
! 157: *> ABSTOL is DOUBLE PRECISION
! 158: *> The absolute error tolerance for the eigenvalues.
! 159: *> An approximate eigenvalue is accepted as converged
! 160: *> when it is determined to lie in an interval [a,b]
! 161: *> of width less than or equal to
! 162: *>
! 163: *> ABSTOL + EPS * max( |a|,|b| ) ,
! 164: *>
! 165: *> where EPS is the machine precision. If ABSTOL is less than
! 166: *> or equal to zero, then EPS*|T| will be used in its place,
! 167: *> where |T| is the 1-norm of the tridiagonal matrix obtained
! 168: *> by reducing A to tridiagonal form.
! 169: *>
! 170: *> See "Computing Small Singular Values of Bidiagonal Matrices
! 171: *> with Guaranteed High Relative Accuracy," by Demmel and
! 172: *> Kahan, LAPACK Working Note #3.
! 173: *>
! 174: *> If high relative accuracy is important, set ABSTOL to
! 175: *> DLAMCH( 'Safe minimum' ). Doing so will guarantee that
! 176: *> eigenvalues are computed to high relative accuracy when
! 177: *> possible in future releases. The current code does not
! 178: *> make any guarantees about high relative accuracy, but
! 179: *> future releases will. See J. Barlow and J. Demmel,
! 180: *> "Computing Accurate Eigensystems of Scaled Diagonally
! 181: *> Dominant Matrices", LAPACK Working Note #7, for a discussion
! 182: *> of which matrices define their eigenvalues to high relative
! 183: *> accuracy.
! 184: *> \endverbatim
! 185: *>
! 186: *> \param[out] M
! 187: *> \verbatim
! 188: *> M is INTEGER
! 189: *> The total number of eigenvalues found. 0 <= M <= N.
! 190: *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
! 191: *> \endverbatim
! 192: *>
! 193: *> \param[out] W
! 194: *> \verbatim
! 195: *> W is DOUBLE PRECISION array, dimension (N)
! 196: *> The first M elements contain the selected eigenvalues in
! 197: *> ascending order.
! 198: *> \endverbatim
! 199: *>
! 200: *> \param[out] Z
! 201: *> \verbatim
! 202: *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
! 203: *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
! 204: *> contain the orthonormal eigenvectors of the matrix A
! 205: *> corresponding to the selected eigenvalues, with the i-th
! 206: *> column of Z holding the eigenvector associated with W(i).
! 207: *> Note: the user must ensure that at least max(1,M) columns are
! 208: *> supplied in the array Z; if RANGE = 'V', the exact value of M
! 209: *> is not known in advance and an upper bound must be used.
! 210: *> \endverbatim
! 211: *>
! 212: *> \param[in] LDZ
! 213: *> \verbatim
! 214: *> LDZ is INTEGER
! 215: *> The leading dimension of the array Z. LDZ >= 1, and if
! 216: *> JOBZ = 'V', LDZ >= max(1,N).
! 217: *> \endverbatim
! 218: *>
! 219: *> \param[out] ISUPPZ
! 220: *> \verbatim
! 221: *> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
! 222: *> The support of the eigenvectors in Z, i.e., the indices
! 223: *> indicating the nonzero elements in Z. The i-th eigenvector
! 224: *> is nonzero only in elements ISUPPZ( 2*i-1 ) through
! 225: *> ISUPPZ( 2*i ).
! 226: *> Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
! 227: *> \endverbatim
! 228: *>
! 229: *> \param[out] WORK
! 230: *> \verbatim
! 231: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 232: *> On exit, if INFO = 0, WORK(1) returns the optimal (and
! 233: *> minimal) LWORK.
! 234: *> \endverbatim
! 235: *>
! 236: *> \param[in] LWORK
! 237: *> \verbatim
! 238: *> LWORK is INTEGER
! 239: *> The dimension of the array WORK. LWORK >= max(1,20*N).
! 240: *>
! 241: *> If LWORK = -1, then a workspace query is assumed; the routine
! 242: *> only calculates the optimal sizes of the WORK and IWORK
! 243: *> arrays, returns these values as the first entries of the WORK
! 244: *> and IWORK arrays, and no error message related to LWORK or
! 245: *> LIWORK is issued by XERBLA.
! 246: *> \endverbatim
! 247: *>
! 248: *> \param[out] IWORK
! 249: *> \verbatim
! 250: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
! 251: *> On exit, if INFO = 0, IWORK(1) returns the optimal (and
! 252: *> minimal) LIWORK.
! 253: *> \endverbatim
! 254: *>
! 255: *> \param[in] LIWORK
! 256: *> \verbatim
! 257: *> LIWORK is INTEGER
! 258: *> The dimension of the array IWORK. LIWORK >= max(1,10*N).
! 259: *>
! 260: *> If LIWORK = -1, then a workspace query is assumed; the
! 261: *> routine only calculates the optimal sizes of the WORK and
! 262: *> IWORK arrays, returns these values as the first entries of
! 263: *> the WORK and IWORK arrays, and no error message related to
! 264: *> LWORK or LIWORK is issued by XERBLA.
! 265: *> \endverbatim
! 266: *>
! 267: *> \param[out] INFO
! 268: *> \verbatim
! 269: *> INFO is INTEGER
! 270: *> = 0: successful exit
! 271: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 272: *> > 0: Internal error
! 273: *> \endverbatim
! 274: *
! 275: * Authors:
! 276: * ========
! 277: *
! 278: *> \author Univ. of Tennessee
! 279: *> \author Univ. of California Berkeley
! 280: *> \author Univ. of Colorado Denver
! 281: *> \author NAG Ltd.
! 282: *
! 283: *> \date November 2011
! 284: *
! 285: *> \ingroup doubleOTHEReigen
! 286: *
! 287: *> \par Contributors:
! 288: * ==================
! 289: *>
! 290: *> Inderjit Dhillon, IBM Almaden, USA \n
! 291: *> Osni Marques, LBNL/NERSC, USA \n
! 292: *> Ken Stanley, Computer Science Division, University of
! 293: *> California at Berkeley, USA \n
! 294: *>
! 295: * =====================================================================
1.1 bertrand 296: SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
297: $ M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK,
298: $ LIWORK, INFO )
299: *
1.8 ! bertrand 300: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 301: * -- LAPACK is a software package provided by Univ. of Tennessee, --
302: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 303: * November 2011
1.1 bertrand 304: *
305: * .. Scalar Arguments ..
306: CHARACTER JOBZ, RANGE
307: INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
308: DOUBLE PRECISION ABSTOL, VL, VU
309: * ..
310: * .. Array Arguments ..
311: INTEGER ISUPPZ( * ), IWORK( * )
312: DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
313: * ..
314: *
315: * =====================================================================
316: *
317: * .. Parameters ..
318: DOUBLE PRECISION ZERO, ONE, TWO
319: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
320: * ..
321: * .. Local Scalars ..
322: LOGICAL ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ,
323: $ TRYRAC
324: CHARACTER ORDER
325: INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP,
326: $ INDIWO, ISCALE, ITMP1, J, JJ, LIWMIN, LWMIN,
327: $ NSPLIT
328: DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
329: $ TMP1, TNRM, VLL, VUU
330: * ..
331: * .. External Functions ..
332: LOGICAL LSAME
333: INTEGER ILAENV
334: DOUBLE PRECISION DLAMCH, DLANST
335: EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
336: * ..
337: * .. External Subroutines ..
338: EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTEMR, DSTEIN, DSTERF,
339: $ DSWAP, XERBLA
340: * ..
341: * .. Intrinsic Functions ..
342: INTRINSIC MAX, MIN, SQRT
343: * ..
344: * .. Executable Statements ..
345: *
346: *
347: * Test the input parameters.
348: *
349: IEEEOK = ILAENV( 10, 'DSTEVR', 'N', 1, 2, 3, 4 )
350: *
351: WANTZ = LSAME( JOBZ, 'V' )
352: ALLEIG = LSAME( RANGE, 'A' )
353: VALEIG = LSAME( RANGE, 'V' )
354: INDEIG = LSAME( RANGE, 'I' )
355: *
356: LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) )
357: LWMIN = MAX( 1, 20*N )
358: LIWMIN = MAX( 1, 10*N )
359: *
360: *
361: INFO = 0
362: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
363: INFO = -1
364: ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
365: INFO = -2
366: ELSE IF( N.LT.0 ) THEN
367: INFO = -3
368: ELSE
369: IF( VALEIG ) THEN
370: IF( N.GT.0 .AND. VU.LE.VL )
371: $ INFO = -7
372: ELSE IF( INDEIG ) THEN
373: IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
374: INFO = -8
375: ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
376: INFO = -9
377: END IF
378: END IF
379: END IF
380: IF( INFO.EQ.0 ) THEN
381: IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
382: INFO = -14
383: END IF
384: END IF
385: *
386: IF( INFO.EQ.0 ) THEN
387: WORK( 1 ) = LWMIN
388: IWORK( 1 ) = LIWMIN
389: *
390: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
391: INFO = -17
392: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
393: INFO = -19
394: END IF
395: END IF
396: *
397: IF( INFO.NE.0 ) THEN
398: CALL XERBLA( 'DSTEVR', -INFO )
399: RETURN
400: ELSE IF( LQUERY ) THEN
401: RETURN
402: END IF
403: *
404: * Quick return if possible
405: *
406: M = 0
407: IF( N.EQ.0 )
408: $ RETURN
409: *
410: IF( N.EQ.1 ) THEN
411: IF( ALLEIG .OR. INDEIG ) THEN
412: M = 1
413: W( 1 ) = D( 1 )
414: ELSE
415: IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN
416: M = 1
417: W( 1 ) = D( 1 )
418: END IF
419: END IF
420: IF( WANTZ )
421: $ Z( 1, 1 ) = ONE
422: RETURN
423: END IF
424: *
425: * Get machine constants.
426: *
427: SAFMIN = DLAMCH( 'Safe minimum' )
428: EPS = DLAMCH( 'Precision' )
429: SMLNUM = SAFMIN / EPS
430: BIGNUM = ONE / SMLNUM
431: RMIN = SQRT( SMLNUM )
432: RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
433: *
434: *
435: * Scale matrix to allowable range, if necessary.
436: *
437: ISCALE = 0
438: VLL = VL
439: VUU = VU
440: *
441: TNRM = DLANST( 'M', N, D, E )
442: IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
443: ISCALE = 1
444: SIGMA = RMIN / TNRM
445: ELSE IF( TNRM.GT.RMAX ) THEN
446: ISCALE = 1
447: SIGMA = RMAX / TNRM
448: END IF
449: IF( ISCALE.EQ.1 ) THEN
450: CALL DSCAL( N, SIGMA, D, 1 )
451: CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
452: IF( VALEIG ) THEN
453: VLL = VL*SIGMA
454: VUU = VU*SIGMA
455: END IF
456: END IF
457:
458: * Initialize indices into workspaces. Note: These indices are used only
459: * if DSTERF or DSTEMR fail.
460:
461: * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
462: * stores the block indices of each of the M<=N eigenvalues.
463: INDIBL = 1
464: * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
465: * stores the starting and finishing indices of each block.
466: INDISP = INDIBL + N
467: * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
468: * that corresponding to eigenvectors that fail to converge in
469: * DSTEIN. This information is discarded; if any fail, the driver
470: * returns INFO > 0.
471: INDIFL = INDISP + N
472: * INDIWO is the offset of the remaining integer workspace.
473: INDIWO = INDISP + N
474: *
475: * If all eigenvalues are desired, then
476: * call DSTERF or DSTEMR. If this fails for some eigenvalue, then
477: * try DSTEBZ.
478: *
479: *
480: TEST = .FALSE.
481: IF( INDEIG ) THEN
482: IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
483: TEST = .TRUE.
484: END IF
485: END IF
486: IF( ( ALLEIG .OR. TEST ) .AND. IEEEOK.EQ.1 ) THEN
487: CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 )
488: IF( .NOT.WANTZ ) THEN
489: CALL DCOPY( N, D, 1, W, 1 )
490: CALL DSTERF( N, W, WORK, INFO )
491: ELSE
492: CALL DCOPY( N, D, 1, WORK( N+1 ), 1 )
493: IF (ABSTOL .LE. TWO*N*EPS) THEN
494: TRYRAC = .TRUE.
495: ELSE
496: TRYRAC = .FALSE.
497: END IF
498: CALL DSTEMR( JOBZ, 'A', N, WORK( N+1 ), WORK, VL, VU, IL,
499: $ IU, M, W, Z, LDZ, N, ISUPPZ, TRYRAC,
500: $ WORK( 2*N+1 ), LWORK-2*N, IWORK, LIWORK, INFO )
501: *
502: END IF
503: IF( INFO.EQ.0 ) THEN
504: M = N
505: GO TO 10
506: END IF
507: INFO = 0
508: END IF
509: *
510: * Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN.
511: *
512: IF( WANTZ ) THEN
513: ORDER = 'B'
514: ELSE
515: ORDER = 'E'
516: END IF
517:
518: CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M,
519: $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), WORK,
520: $ IWORK( INDIWO ), INFO )
521: *
522: IF( WANTZ ) THEN
523: CALL DSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ),
524: $ Z, LDZ, WORK, IWORK( INDIWO ), IWORK( INDIFL ),
525: $ INFO )
526: END IF
527: *
528: * If matrix was scaled, then rescale eigenvalues appropriately.
529: *
530: 10 CONTINUE
531: IF( ISCALE.EQ.1 ) THEN
532: IF( INFO.EQ.0 ) THEN
533: IMAX = M
534: ELSE
535: IMAX = INFO - 1
536: END IF
537: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
538: END IF
539: *
540: * If eigenvalues are not in order, then sort them, along with
541: * eigenvectors.
542: *
543: IF( WANTZ ) THEN
544: DO 30 J = 1, M - 1
545: I = 0
546: TMP1 = W( J )
547: DO 20 JJ = J + 1, M
548: IF( W( JJ ).LT.TMP1 ) THEN
549: I = JJ
550: TMP1 = W( JJ )
551: END IF
552: 20 CONTINUE
553: *
554: IF( I.NE.0 ) THEN
555: ITMP1 = IWORK( I )
556: W( I ) = W( J )
557: IWORK( I ) = IWORK( J )
558: W( J ) = TMP1
559: IWORK( J ) = ITMP1
560: CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
561: END IF
562: 30 CONTINUE
563: END IF
564: *
565: * Causes problems with tests 19 & 20:
566: * IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002
567: *
568: *
569: WORK( 1 ) = LWMIN
570: IWORK( 1 ) = LIWMIN
571: RETURN
572: *
573: * End of DSTEVR
574: *
575: END
CVSweb interface <joel.bertrand@systella.fr>