Annotation of rpl/lapack/lapack/dsbevx_2stage.f, revision 1.1

1.1     ! bertrand    1: *> \brief <b> DSBEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
        !             2: *
        !             3: *  @precisions fortran d -> s
        !             4: *
        !             5: *  =========== DOCUMENTATION ===========
        !             6: *
        !             7: * Online html documentation available at
        !             8: *            http://www.netlib.org/lapack/explore-html/
        !             9: *
        !            10: *> \htmlonly
        !            11: *> Download DSBEVX_2STAGE + dependencies
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbevx_2stage.f">
        !            13: *> [TGZ]</a>
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbevx_2stage.f">
        !            15: *> [ZIP]</a>
        !            16: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbevx_2stage.f">
        !            17: *> [TXT]</a>
        !            18: *> \endhtmlonly
        !            19: *
        !            20: *  Definition:
        !            21: *  ===========
        !            22: *
        !            23: *       SUBROUTINE DSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q,
        !            24: *                                 LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z,
        !            25: *                                 LDZ, WORK, LWORK, IWORK, IFAIL, INFO )
        !            26: *
        !            27: *       IMPLICIT NONE
        !            28: *
        !            29: *       .. Scalar Arguments ..
        !            30: *       CHARACTER          JOBZ, RANGE, UPLO
        !            31: *       INTEGER            IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N, LWORK
        !            32: *       DOUBLE PRECISION   ABSTOL, VL, VU
        !            33: *       ..
        !            34: *       .. Array Arguments ..
        !            35: *       INTEGER            IFAIL( * ), IWORK( * )
        !            36: *       DOUBLE PRECISION   AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
        !            37: *      $                   Z( LDZ, * )
        !            38: *       ..
        !            39: *
        !            40: *
        !            41: *> \par Purpose:
        !            42: *  =============
        !            43: *>
        !            44: *> \verbatim
        !            45: *>
        !            46: *> DSBEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
        !            47: *> of a real symmetric band matrix A using the 2stage technique for
        !            48: *> the reduction to tridiagonal. Eigenvalues and eigenvectors can
        !            49: *> be selected by specifying either a range of values or a range of
        !            50: *> indices for the desired eigenvalues.
        !            51: *> \endverbatim
        !            52: *
        !            53: *  Arguments:
        !            54: *  ==========
        !            55: *
        !            56: *> \param[in] JOBZ
        !            57: *> \verbatim
        !            58: *>          JOBZ is CHARACTER*1
        !            59: *>          = 'N':  Compute eigenvalues only;
        !            60: *>          = 'V':  Compute eigenvalues and eigenvectors.
        !            61: *>                  Not available in this release.
        !            62: *> \endverbatim
        !            63: *>
        !            64: *> \param[in] RANGE
        !            65: *> \verbatim
        !            66: *>          RANGE is CHARACTER*1
        !            67: *>          = 'A': all eigenvalues will be found;
        !            68: *>          = 'V': all eigenvalues in the half-open interval (VL,VU]
        !            69: *>                 will be found;
        !            70: *>          = 'I': the IL-th through IU-th eigenvalues will be found.
        !            71: *> \endverbatim
        !            72: *>
        !            73: *> \param[in] UPLO
        !            74: *> \verbatim
        !            75: *>          UPLO is CHARACTER*1
        !            76: *>          = 'U':  Upper triangle of A is stored;
        !            77: *>          = 'L':  Lower triangle of A is stored.
        !            78: *> \endverbatim
        !            79: *>
        !            80: *> \param[in] N
        !            81: *> \verbatim
        !            82: *>          N is INTEGER
        !            83: *>          The order of the matrix A.  N >= 0.
        !            84: *> \endverbatim
        !            85: *>
        !            86: *> \param[in] KD
        !            87: *> \verbatim
        !            88: *>          KD is INTEGER
        !            89: *>          The number of superdiagonals of the matrix A if UPLO = 'U',
        !            90: *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
        !            91: *> \endverbatim
        !            92: *>
        !            93: *> \param[in,out] AB
        !            94: *> \verbatim
        !            95: *>          AB is DOUBLE PRECISION array, dimension (LDAB, N)
        !            96: *>          On entry, the upper or lower triangle of the symmetric band
        !            97: *>          matrix A, stored in the first KD+1 rows of the array.  The
        !            98: *>          j-th column of A is stored in the j-th column of the array AB
        !            99: *>          as follows:
        !           100: *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
        !           101: *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
        !           102: *>
        !           103: *>          On exit, AB is overwritten by values generated during the
        !           104: *>          reduction to tridiagonal form.  If UPLO = 'U', the first
        !           105: *>          superdiagonal and the diagonal of the tridiagonal matrix T
        !           106: *>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
        !           107: *>          the diagonal and first subdiagonal of T are returned in the
        !           108: *>          first two rows of AB.
        !           109: *> \endverbatim
        !           110: *>
        !           111: *> \param[in] LDAB
        !           112: *> \verbatim
        !           113: *>          LDAB is INTEGER
        !           114: *>          The leading dimension of the array AB.  LDAB >= KD + 1.
        !           115: *> \endverbatim
        !           116: *>
        !           117: *> \param[out] Q
        !           118: *> \verbatim
        !           119: *>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
        !           120: *>          If JOBZ = 'V', the N-by-N orthogonal matrix used in the
        !           121: *>                         reduction to tridiagonal form.
        !           122: *>          If JOBZ = 'N', the array Q is not referenced.
        !           123: *> \endverbatim
        !           124: *>
        !           125: *> \param[in] LDQ
        !           126: *> \verbatim
        !           127: *>          LDQ is INTEGER
        !           128: *>          The leading dimension of the array Q.  If JOBZ = 'V', then
        !           129: *>          LDQ >= max(1,N).
        !           130: *> \endverbatim
        !           131: *>
        !           132: *> \param[in] VL
        !           133: *> \verbatim
        !           134: *>          VL is DOUBLE PRECISION
        !           135: *>          If RANGE='V', the lower bound of the interval to
        !           136: *>          be searched for eigenvalues. VL < VU.
        !           137: *>          Not referenced if RANGE = 'A' or 'I'.
        !           138: *> \endverbatim
        !           139: *>
        !           140: *> \param[in] VU
        !           141: *> \verbatim
        !           142: *>          VU is DOUBLE PRECISION
        !           143: *>          If RANGE='V', the upper bound of the interval to
        !           144: *>          be searched for eigenvalues. VL < VU.
        !           145: *>          Not referenced if RANGE = 'A' or 'I'.
        !           146: *> \endverbatim
        !           147: *>
        !           148: *> \param[in] IL
        !           149: *> \verbatim
        !           150: *>          IL is INTEGER
        !           151: *>          If RANGE='I', the index of the
        !           152: *>          smallest eigenvalue to be returned.
        !           153: *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
        !           154: *>          Not referenced if RANGE = 'A' or 'V'.
        !           155: *> \endverbatim
        !           156: *>
        !           157: *> \param[in] IU
        !           158: *> \verbatim
        !           159: *>          IU is INTEGER
        !           160: *>          If RANGE='I', the index of the
        !           161: *>          largest eigenvalue to be returned.
        !           162: *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
        !           163: *>          Not referenced if RANGE = 'A' or 'V'.
        !           164: *> \endverbatim
        !           165: *>
        !           166: *> \param[in] ABSTOL
        !           167: *> \verbatim
        !           168: *>          ABSTOL is DOUBLE PRECISION
        !           169: *>          The absolute error tolerance for the eigenvalues.
        !           170: *>          An approximate eigenvalue is accepted as converged
        !           171: *>          when it is determined to lie in an interval [a,b]
        !           172: *>          of width less than or equal to
        !           173: *>
        !           174: *>                  ABSTOL + EPS *   max( |a|,|b| ) ,
        !           175: *>
        !           176: *>          where EPS is the machine precision.  If ABSTOL is less than
        !           177: *>          or equal to zero, then  EPS*|T|  will be used in its place,
        !           178: *>          where |T| is the 1-norm of the tridiagonal matrix obtained
        !           179: *>          by reducing AB to tridiagonal form.
        !           180: *>
        !           181: *>          Eigenvalues will be computed most accurately when ABSTOL is
        !           182: *>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
        !           183: *>          If this routine returns with INFO>0, indicating that some
        !           184: *>          eigenvectors did not converge, try setting ABSTOL to
        !           185: *>          2*DLAMCH('S').
        !           186: *>
        !           187: *>          See "Computing Small Singular Values of Bidiagonal Matrices
        !           188: *>          with Guaranteed High Relative Accuracy," by Demmel and
        !           189: *>          Kahan, LAPACK Working Note #3.
        !           190: *> \endverbatim
        !           191: *>
        !           192: *> \param[out] M
        !           193: *> \verbatim
        !           194: *>          M is INTEGER
        !           195: *>          The total number of eigenvalues found.  0 <= M <= N.
        !           196: *>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
        !           197: *> \endverbatim
        !           198: *>
        !           199: *> \param[out] W
        !           200: *> \verbatim
        !           201: *>          W is DOUBLE PRECISION array, dimension (N)
        !           202: *>          The first M elements contain the selected eigenvalues in
        !           203: *>          ascending order.
        !           204: *> \endverbatim
        !           205: *>
        !           206: *> \param[out] Z
        !           207: *> \verbatim
        !           208: *>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
        !           209: *>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
        !           210: *>          contain the orthonormal eigenvectors of the matrix A
        !           211: *>          corresponding to the selected eigenvalues, with the i-th
        !           212: *>          column of Z holding the eigenvector associated with W(i).
        !           213: *>          If an eigenvector fails to converge, then that column of Z
        !           214: *>          contains the latest approximation to the eigenvector, and the
        !           215: *>          index of the eigenvector is returned in IFAIL.
        !           216: *>          If JOBZ = 'N', then Z is not referenced.
        !           217: *>          Note: the user must ensure that at least max(1,M) columns are
        !           218: *>          supplied in the array Z; if RANGE = 'V', the exact value of M
        !           219: *>          is not known in advance and an upper bound must be used.
        !           220: *> \endverbatim
        !           221: *>
        !           222: *> \param[in] LDZ
        !           223: *> \verbatim
        !           224: *>          LDZ is INTEGER
        !           225: *>          The leading dimension of the array Z.  LDZ >= 1, and if
        !           226: *>          JOBZ = 'V', LDZ >= max(1,N).
        !           227: *> \endverbatim
        !           228: *>
        !           229: *> \param[out] WORK
        !           230: *> \verbatim
        !           231: *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
        !           232: *> \endverbatim
        !           233: *>
        !           234: *> \param[in] LWORK
        !           235: *> \verbatim
        !           236: *>          LWORK is INTEGER
        !           237: *>          The length of the array WORK. LWORK >= 1, when N <= 1;
        !           238: *>          otherwise  
        !           239: *>          If JOBZ = 'N' and N > 1, LWORK must be queried.
        !           240: *>                                   LWORK = MAX(1, 7*N, dimension) where
        !           241: *>                                   dimension = (2KD+1)*N + KD*NTHREADS + 2*N
        !           242: *>                                   where KD is the size of the band.
        !           243: *>                                   NTHREADS is the number of threads used when
        !           244: *>                                   openMP compilation is enabled, otherwise =1.
        !           245: *>          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available
        !           246: *>
        !           247: *>          If LWORK = -1, then a workspace query is assumed; the routine
        !           248: *>          only calculates the optimal size of the WORK array, returns
        !           249: *>          this value as the first entry of the WORK array, and no error
        !           250: *>          message related to LWORK is issued by XERBLA.
        !           251: *> \endverbatim
        !           252: *>
        !           253: *> \param[out] IWORK
        !           254: *> \verbatim
        !           255: *>          IWORK is INTEGER array, dimension (5*N)
        !           256: *> \endverbatim
        !           257: *>
        !           258: *> \param[out] IFAIL
        !           259: *> \verbatim
        !           260: *>          IFAIL is INTEGER array, dimension (N)
        !           261: *>          If JOBZ = 'V', then if INFO = 0, the first M elements of
        !           262: *>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
        !           263: *>          indices of the eigenvectors that failed to converge.
        !           264: *>          If JOBZ = 'N', then IFAIL is not referenced.
        !           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:  if INFO = i, then i eigenvectors failed to converge.
        !           273: *>                Their indices are stored in array IFAIL.
        !           274: *> \endverbatim
        !           275: *
        !           276: *  Authors:
        !           277: *  ========
        !           278: *
        !           279: *> \author Univ. of Tennessee
        !           280: *> \author Univ. of California Berkeley
        !           281: *> \author Univ. of Colorado Denver
        !           282: *> \author NAG Ltd.
        !           283: *
        !           284: *> \date June 2016
        !           285: *
        !           286: *> \ingroup doubleOTHEReigen
        !           287: *
        !           288: *> \par Further Details:
        !           289: *  =====================
        !           290: *>
        !           291: *> \verbatim
        !           292: *>
        !           293: *>  All details about the 2stage techniques are available in:
        !           294: *>
        !           295: *>  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
        !           296: *>  Parallel reduction to condensed forms for symmetric eigenvalue problems
        !           297: *>  using aggregated fine-grained and memory-aware kernels. In Proceedings
        !           298: *>  of 2011 International Conference for High Performance Computing,
        !           299: *>  Networking, Storage and Analysis (SC '11), New York, NY, USA,
        !           300: *>  Article 8 , 11 pages.
        !           301: *>  http://doi.acm.org/10.1145/2063384.2063394
        !           302: *>
        !           303: *>  A. Haidar, J. Kurzak, P. Luszczek, 2013.
        !           304: *>  An improved parallel singular value algorithm and its implementation 
        !           305: *>  for multicore hardware, In Proceedings of 2013 International Conference
        !           306: *>  for High Performance Computing, Networking, Storage and Analysis (SC '13).
        !           307: *>  Denver, Colorado, USA, 2013.
        !           308: *>  Article 90, 12 pages.
        !           309: *>  http://doi.acm.org/10.1145/2503210.2503292
        !           310: *>
        !           311: *>  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
        !           312: *>  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
        !           313: *>  calculations based on fine-grained memory aware tasks.
        !           314: *>  International Journal of High Performance Computing Applications.
        !           315: *>  Volume 28 Issue 2, Pages 196-209, May 2014.
        !           316: *>  http://hpc.sagepub.com/content/28/2/196 
        !           317: *>
        !           318: *> \endverbatim
        !           319: *
        !           320: *  =====================================================================
        !           321:       SUBROUTINE DSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q,
        !           322:      $                          LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z,
        !           323:      $                          LDZ, WORK, LWORK, IWORK, IFAIL, INFO )
        !           324: *
        !           325:       IMPLICIT NONE
        !           326: *
        !           327: *  -- LAPACK driver routine (version 3.7.0) --
        !           328: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
        !           329: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
        !           330: *     June 2016
        !           331: *
        !           332: *     .. Scalar Arguments ..
        !           333:       CHARACTER          JOBZ, RANGE, UPLO
        !           334:       INTEGER            IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N, LWORK
        !           335:       DOUBLE PRECISION   ABSTOL, VL, VU
        !           336: *     ..
        !           337: *     .. Array Arguments ..
        !           338:       INTEGER            IFAIL( * ), IWORK( * )
        !           339:       DOUBLE PRECISION   AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
        !           340:      $                   Z( LDZ, * )
        !           341: *     ..
        !           342: *
        !           343: *  =====================================================================
        !           344: *
        !           345: *     .. Parameters ..
        !           346:       DOUBLE PRECISION   ZERO, ONE
        !           347:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
        !           348: *     ..
        !           349: *     .. Local Scalars ..
        !           350:       LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ,
        !           351:      $                   LQUERY
        !           352:       CHARACTER          ORDER
        !           353:       INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
        !           354:      $                   INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
        !           355:      $                   LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
        !           356:      $                   NSPLIT
        !           357:       DOUBLE PRECISION   ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
        !           358:      $                   SIGMA, SMLNUM, TMP1, VLL, VUU
        !           359: *     ..
        !           360: *     .. External Functions ..
        !           361:       LOGICAL            LSAME
        !           362:       INTEGER            ILAENV
        !           363:       DOUBLE PRECISION   DLAMCH, DLANSB
        !           364:       EXTERNAL           LSAME, DLAMCH, DLANSB, ILAENV
        !           365: *     ..
        !           366: *     .. External Subroutines ..
        !           367:       EXTERNAL           DCOPY, DGEMV, DLACPY, DLASCL, DSCAL,
        !           368:      $                   DSTEBZ, DSTEIN, DSTEQR, DSTERF, DSWAP, XERBLA,
        !           369:      $                   DSYTRD_SB2ST
        !           370: *     ..
        !           371: *     .. Intrinsic Functions ..
        !           372:       INTRINSIC          MAX, MIN, SQRT
        !           373: *     ..
        !           374: *     .. Executable Statements ..
        !           375: *
        !           376: *     Test the input parameters.
        !           377: *
        !           378:       WANTZ = LSAME( JOBZ, 'V' )
        !           379:       ALLEIG = LSAME( RANGE, 'A' )
        !           380:       VALEIG = LSAME( RANGE, 'V' )
        !           381:       INDEIG = LSAME( RANGE, 'I' )
        !           382:       LOWER = LSAME( UPLO, 'L' )
        !           383:       LQUERY = ( LWORK.EQ.-1 )
        !           384: *
        !           385:       INFO = 0
        !           386:       IF( .NOT.( LSAME( JOBZ, 'N' ) ) ) THEN
        !           387:          INFO = -1
        !           388:       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
        !           389:          INFO = -2
        !           390:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
        !           391:          INFO = -3
        !           392:       ELSE IF( N.LT.0 ) THEN
        !           393:          INFO = -4
        !           394:       ELSE IF( KD.LT.0 ) THEN
        !           395:          INFO = -5
        !           396:       ELSE IF( LDAB.LT.KD+1 ) THEN
        !           397:          INFO = -7
        !           398:       ELSE IF( WANTZ .AND. LDQ.LT.MAX( 1, N ) ) THEN
        !           399:          INFO = -9
        !           400:       ELSE
        !           401:          IF( VALEIG ) THEN
        !           402:             IF( N.GT.0 .AND. VU.LE.VL )
        !           403:      $         INFO = -11
        !           404:          ELSE IF( INDEIG ) THEN
        !           405:             IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
        !           406:                INFO = -12
        !           407:             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
        !           408:                INFO = -13
        !           409:             END IF
        !           410:          END IF
        !           411:       END IF
        !           412:       IF( INFO.EQ.0 ) THEN
        !           413:          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) )
        !           414:      $      INFO = -18
        !           415:       END IF
        !           416: *
        !           417:       IF( INFO.EQ.0 ) THEN
        !           418:          IF( N.LE.1 ) THEN
        !           419:             LWMIN = 1
        !           420:             WORK( 1 ) = LWMIN
        !           421:          ELSE
        !           422:             IB    = ILAENV( 18, 'DSYTRD_SB2ST', JOBZ, N, KD, -1, -1 )
        !           423:             LHTRD = ILAENV( 19, 'DSYTRD_SB2ST', JOBZ, N, KD, IB, -1 )
        !           424:             LWTRD = ILAENV( 20, 'DSYTRD_SB2ST', JOBZ, N, KD, IB, -1 )
        !           425:             LWMIN = 2*N + LHTRD + LWTRD
        !           426:             WORK( 1 )  = LWMIN
        !           427:          ENDIF
        !           428: *
        !           429:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY )
        !           430:      $      INFO = -20
        !           431:       END IF
        !           432: *
        !           433:       IF( INFO.NE.0 ) THEN
        !           434:          CALL XERBLA( 'DSBEVX_2STAGE ', -INFO )
        !           435:          RETURN
        !           436:       ELSE IF( LQUERY ) THEN
        !           437:          RETURN
        !           438:       END IF
        !           439: *
        !           440: *     Quick return if possible
        !           441: *
        !           442:       M = 0
        !           443:       IF( N.EQ.0 )
        !           444:      $   RETURN
        !           445: *
        !           446:       IF( N.EQ.1 ) THEN
        !           447:          M = 1
        !           448:          IF( LOWER ) THEN
        !           449:             TMP1 = AB( 1, 1 )
        !           450:          ELSE
        !           451:             TMP1 = AB( KD+1, 1 )
        !           452:          END IF
        !           453:          IF( VALEIG ) THEN
        !           454:             IF( .NOT.( VL.LT.TMP1 .AND. VU.GE.TMP1 ) )
        !           455:      $         M = 0
        !           456:          END IF
        !           457:          IF( M.EQ.1 ) THEN
        !           458:             W( 1 ) = TMP1
        !           459:             IF( WANTZ )
        !           460:      $         Z( 1, 1 ) = ONE
        !           461:          END IF
        !           462:          RETURN
        !           463:       END IF
        !           464: *
        !           465: *     Get machine constants.
        !           466: *
        !           467:       SAFMIN = DLAMCH( 'Safe minimum' )
        !           468:       EPS    = DLAMCH( 'Precision' )
        !           469:       SMLNUM = SAFMIN / EPS
        !           470:       BIGNUM = ONE / SMLNUM
        !           471:       RMIN   = SQRT( SMLNUM )
        !           472:       RMAX   = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
        !           473: *
        !           474: *     Scale matrix to allowable range, if necessary.
        !           475: *
        !           476:       ISCALE = 0
        !           477:       ABSTLL = ABSTOL
        !           478:       IF( VALEIG ) THEN
        !           479:          VLL = VL
        !           480:          VUU = VU
        !           481:       ELSE
        !           482:          VLL = ZERO
        !           483:          VUU = ZERO
        !           484:       END IF
        !           485:       ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
        !           486:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
        !           487:          ISCALE = 1
        !           488:          SIGMA = RMIN / ANRM
        !           489:       ELSE IF( ANRM.GT.RMAX ) THEN
        !           490:          ISCALE = 1
        !           491:          SIGMA = RMAX / ANRM
        !           492:       END IF
        !           493:       IF( ISCALE.EQ.1 ) THEN
        !           494:          IF( LOWER ) THEN
        !           495:             CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
        !           496:          ELSE
        !           497:             CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
        !           498:          END IF
        !           499:          IF( ABSTOL.GT.0 )
        !           500:      $      ABSTLL = ABSTOL*SIGMA
        !           501:          IF( VALEIG ) THEN
        !           502:             VLL = VL*SIGMA
        !           503:             VUU = VU*SIGMA
        !           504:          END IF
        !           505:       END IF
        !           506: *
        !           507: *     Call DSYTRD_SB2ST to reduce symmetric band matrix to tridiagonal form.
        !           508: *
        !           509:       INDD    = 1
        !           510:       INDE    = INDD + N
        !           511:       INDHOUS = INDE + N
        !           512:       INDWRK  = INDHOUS + LHTRD
        !           513:       LLWORK  = LWORK - INDWRK + 1
        !           514: *
        !           515:       CALL DSYTRD_SB2ST( "N", JOBZ, UPLO, N, KD, AB, LDAB, WORK( INDD ),
        !           516:      $                    WORK( INDE ), WORK( INDHOUS ), LHTRD, 
        !           517:      $                    WORK( INDWRK ), LLWORK, IINFO )
        !           518: *
        !           519: *     If all eigenvalues are desired and ABSTOL is less than or equal
        !           520: *     to zero, then call DSTERF or SSTEQR.  If this fails for some
        !           521: *     eigenvalue, then try DSTEBZ.
        !           522: *
        !           523:       TEST = .FALSE.
        !           524:       IF (INDEIG) THEN
        !           525:          IF (IL.EQ.1 .AND. IU.EQ.N) THEN
        !           526:             TEST = .TRUE.
        !           527:          END IF
        !           528:       END IF
        !           529:       IF ((ALLEIG .OR. TEST) .AND. (ABSTOL.LE.ZERO)) THEN
        !           530:          CALL DCOPY( N, WORK( INDD ), 1, W, 1 )
        !           531:          INDEE = INDWRK + 2*N
        !           532:          IF( .NOT.WANTZ ) THEN
        !           533:             CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
        !           534:             CALL DSTERF( N, W, WORK( INDEE ), INFO )
        !           535:          ELSE
        !           536:             CALL DLACPY( 'A', N, N, Q, LDQ, Z, LDZ )
        !           537:             CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
        !           538:             CALL DSTEQR( JOBZ, N, W, WORK( INDEE ), Z, LDZ,
        !           539:      $                   WORK( INDWRK ), INFO )
        !           540:             IF( INFO.EQ.0 ) THEN
        !           541:                DO 10 I = 1, N
        !           542:                   IFAIL( I ) = 0
        !           543:    10          CONTINUE
        !           544:             END IF
        !           545:          END IF
        !           546:          IF( INFO.EQ.0 ) THEN
        !           547:             M = N
        !           548:             GO TO 30
        !           549:          END IF
        !           550:          INFO = 0
        !           551:       END IF
        !           552: *
        !           553: *     Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
        !           554: *
        !           555:       IF( WANTZ ) THEN
        !           556:          ORDER = 'B'
        !           557:       ELSE
        !           558:          ORDER = 'E'
        !           559:       END IF
        !           560:       INDIBL = 1
        !           561:       INDISP = INDIBL + N
        !           562:       INDIWO = INDISP + N
        !           563:       CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
        !           564:      $             WORK( INDD ), WORK( INDE ), M, NSPLIT, W,
        !           565:      $             IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWRK ),
        !           566:      $             IWORK( INDIWO ), INFO )
        !           567: *
        !           568:       IF( WANTZ ) THEN
        !           569:          CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W,
        !           570:      $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
        !           571:      $                WORK( INDWRK ), IWORK( INDIWO ), IFAIL, INFO )
        !           572: *
        !           573: *        Apply orthogonal matrix used in reduction to tridiagonal
        !           574: *        form to eigenvectors returned by DSTEIN.
        !           575: *
        !           576:          DO 20 J = 1, M
        !           577:             CALL DCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 )
        !           578:             CALL DGEMV( 'N', N, N, ONE, Q, LDQ, WORK, 1, ZERO,
        !           579:      $                  Z( 1, J ), 1 )
        !           580:    20    CONTINUE
        !           581:       END IF
        !           582: *
        !           583: *     If matrix was scaled, then rescale eigenvalues appropriately.
        !           584: *
        !           585:    30 CONTINUE
        !           586:       IF( ISCALE.EQ.1 ) THEN
        !           587:          IF( INFO.EQ.0 ) THEN
        !           588:             IMAX = M
        !           589:          ELSE
        !           590:             IMAX = INFO - 1
        !           591:          END IF
        !           592:          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
        !           593:       END IF
        !           594: *
        !           595: *     If eigenvalues are not in order, then sort them, along with
        !           596: *     eigenvectors.
        !           597: *
        !           598:       IF( WANTZ ) THEN
        !           599:          DO 50 J = 1, M - 1
        !           600:             I = 0
        !           601:             TMP1 = W( J )
        !           602:             DO 40 JJ = J + 1, M
        !           603:                IF( W( JJ ).LT.TMP1 ) THEN
        !           604:                   I = JJ
        !           605:                   TMP1 = W( JJ )
        !           606:                END IF
        !           607:    40       CONTINUE
        !           608: *
        !           609:             IF( I.NE.0 ) THEN
        !           610:                ITMP1 = IWORK( INDIBL+I-1 )
        !           611:                W( I ) = W( J )
        !           612:                IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
        !           613:                W( J ) = TMP1
        !           614:                IWORK( INDIBL+J-1 ) = ITMP1
        !           615:                CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
        !           616:                IF( INFO.NE.0 ) THEN
        !           617:                   ITMP1 = IFAIL( I )
        !           618:                   IFAIL( I ) = IFAIL( J )
        !           619:                   IFAIL( J ) = ITMP1
        !           620:                END IF
        !           621:             END IF
        !           622:    50    CONTINUE
        !           623:       END IF
        !           624: *
        !           625: *     Set WORK(1) to optimal workspace size.
        !           626: *
        !           627:       WORK( 1 ) = LWMIN
        !           628: *
        !           629:       RETURN
        !           630: *
        !           631: *     End of DSBEVX_2STAGE
        !           632: *
        !           633:       END

CVSweb interface <joel.bertrand@systella.fr>