File:  [local] / rpl / lapack / lapack / dsyevr_2stage.f
Revision 1.1: download - view: text, annotated - select for diffs - revision graph
Sat Jun 17 11:02:51 2017 UTC (6 years, 10 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Ajout des nouveaux fichiers pour lapack 3.7.0.

    1: *> \brief <b> DSYEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY 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 DSYEVR_2STAGE + dependencies
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyevr_2stage.f">
   13: *> [TGZ]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyevr_2stage.f">
   15: *> [ZIP]</a>
   16: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyevr_2stage.f">
   17: *> [TXT]</a>
   18: *> \endhtmlonly
   19: *
   20: *  Definition:
   21: *  ===========
   22: *
   23: *       SUBROUTINE DSYEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
   24: *                          IL, IU, ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK,
   25: *                          LWORK, IWORK, LIWORK, INFO )
   26: *
   27: *       IMPLICIT NONE
   28: *
   29: *       .. Scalar Arguments ..
   30: *       CHARACTER          JOBZ, RANGE, UPLO
   31: *       INTEGER            IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
   32: *       DOUBLE PRECISION   ABSTOL, VL, VU
   33: *       ..
   34: *       .. Array Arguments ..
   35: *       INTEGER            ISUPPZ( * ), IWORK( * )
   36: *       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
   37: *       ..
   38: *
   39: *
   40: *> \par Purpose:
   41: *  =============
   42: *>
   43: *> \verbatim
   44: *>
   45: *> DSYEVR_2STAGE computes selected eigenvalues and, optionally, eigenvectors
   46: *> of a real symmetric matrix A using the 2stage technique for
   47: *> the reduction to tridiagonal.  Eigenvalues and eigenvectors can be
   48: *> selected by specifying either a range of values or a range of
   49: *> indices for the desired eigenvalues.
   50: *>
   51: *> DSYEVR_2STAGE first reduces the matrix A to tridiagonal form T with a call
   52: *> to DSYTRD.  Then, whenever possible, DSYEVR_2STAGE calls DSTEMR to compute
   53: *> the eigenspectrum using Relatively Robust Representations.  DSTEMR
   54: *> computes eigenvalues by the dqds algorithm, while orthogonal
   55: *> eigenvectors are computed from various "good" L D L^T representations
   56: *> (also known as Relatively Robust Representations). Gram-Schmidt
   57: *> orthogonalization is avoided as far as possible. More specifically,
   58: *> the various steps of the algorithm are as follows.
   59: *>
   60: *> For each unreduced block (submatrix) of T,
   61: *>    (a) Compute T - sigma I  = L D L^T, so that L and D
   62: *>        define all the wanted eigenvalues to high relative accuracy.
   63: *>        This means that small relative changes in the entries of D and L
   64: *>        cause only small relative changes in the eigenvalues and
   65: *>        eigenvectors. The standard (unfactored) representation of the
   66: *>        tridiagonal matrix T does not have this property in general.
   67: *>    (b) Compute the eigenvalues to suitable accuracy.
   68: *>        If the eigenvectors are desired, the algorithm attains full
   69: *>        accuracy of the computed eigenvalues only right before
   70: *>        the corresponding vectors have to be computed, see steps c) and d).
   71: *>    (c) For each cluster of close eigenvalues, select a new
   72: *>        shift close to the cluster, find a new factorization, and refine
   73: *>        the shifted eigenvalues to suitable accuracy.
   74: *>    (d) For each eigenvalue with a large enough relative separation compute
   75: *>        the corresponding eigenvector by forming a rank revealing twisted
   76: *>        factorization. Go back to (c) for any clusters that remain.
   77: *>
   78: *> The desired accuracy of the output can be specified by the input
   79: *> parameter ABSTOL.
   80: *>
   81: *> For more details, see DSTEMR's documentation and:
   82: *> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
   83: *>   to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
   84: *>   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
   85: *> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
   86: *>   Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
   87: *>   2004.  Also LAPACK Working Note 154.
   88: *> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
   89: *>   tridiagonal eigenvalue/eigenvector problem",
   90: *>   Computer Science Division Technical Report No. UCB/CSD-97-971,
   91: *>   UC Berkeley, May 1997.
   92: *>
   93: *>
   94: *> Note 1 : DSYEVR_2STAGE calls DSTEMR when the full spectrum is requested
   95: *> on machines which conform to the ieee-754 floating point standard.
   96: *> DSYEVR_2STAGE calls DSTEBZ and SSTEIN on non-ieee machines and
   97: *> when partial spectrum requests are made.
   98: *>
   99: *> Normal execution of DSTEMR may create NaNs and infinities and
  100: *> hence may abort due to a floating point exception in environments
  101: *> which do not handle NaNs and infinities in the ieee standard default
  102: *> manner.
  103: *> \endverbatim
  104: *
  105: *  Arguments:
  106: *  ==========
  107: *
  108: *> \param[in] JOBZ
  109: *> \verbatim
  110: *>          JOBZ is CHARACTER*1
  111: *>          = 'N':  Compute eigenvalues only;
  112: *>          = 'V':  Compute eigenvalues and eigenvectors.
  113: *>                  Not available in this release.
  114: *> \endverbatim
  115: *>
  116: *> \param[in] RANGE
  117: *> \verbatim
  118: *>          RANGE is CHARACTER*1
  119: *>          = 'A': all eigenvalues will be found.
  120: *>          = 'V': all eigenvalues in the half-open interval (VL,VU]
  121: *>                 will be found.
  122: *>          = 'I': the IL-th through IU-th eigenvalues will be found.
  123: *>          For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
  124: *>          DSTEIN are called
  125: *> \endverbatim
  126: *>
  127: *> \param[in] UPLO
  128: *> \verbatim
  129: *>          UPLO is CHARACTER*1
  130: *>          = 'U':  Upper triangle of A is stored;
  131: *>          = 'L':  Lower triangle of A is stored.
  132: *> \endverbatim
  133: *>
  134: *> \param[in] N
  135: *> \verbatim
  136: *>          N is INTEGER
  137: *>          The order of the matrix A.  N >= 0.
  138: *> \endverbatim
  139: *>
  140: *> \param[in,out] A
  141: *> \verbatim
  142: *>          A is DOUBLE PRECISION array, dimension (LDA, N)
  143: *>          On entry, the symmetric matrix A.  If UPLO = 'U', the
  144: *>          leading N-by-N upper triangular part of A contains the
  145: *>          upper triangular part of the matrix A.  If UPLO = 'L',
  146: *>          the leading N-by-N lower triangular part of A contains
  147: *>          the lower triangular part of the matrix A.
  148: *>          On exit, the lower triangle (if UPLO='L') or the upper
  149: *>          triangle (if UPLO='U') of A, including the diagonal, is
  150: *>          destroyed.
  151: *> \endverbatim
  152: *>
  153: *> \param[in] LDA
  154: *> \verbatim
  155: *>          LDA is INTEGER
  156: *>          The leading dimension of the array A.  LDA >= max(1,N).
  157: *> \endverbatim
  158: *>
  159: *> \param[in] VL
  160: *> \verbatim
  161: *>          VL is DOUBLE PRECISION
  162: *>          If RANGE='V', the lower bound of the interval to
  163: *>          be searched for eigenvalues. VL < VU.
  164: *>          Not referenced if RANGE = 'A' or 'I'.
  165: *> \endverbatim
  166: *>
  167: *> \param[in] VU
  168: *> \verbatim
  169: *>          VU is DOUBLE PRECISION
  170: *>          If RANGE='V', the upper bound of the interval to
  171: *>          be searched for eigenvalues. VL < VU.
  172: *>          Not referenced if RANGE = 'A' or 'I'.
  173: *> \endverbatim
  174: *>
  175: *> \param[in] IL
  176: *> \verbatim
  177: *>          IL is INTEGER
  178: *>          If RANGE='I', the index of the
  179: *>          smallest eigenvalue to be returned.
  180: *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
  181: *>          Not referenced if RANGE = 'A' or 'V'.
  182: *> \endverbatim
  183: *>
  184: *> \param[in] IU
  185: *> \verbatim
  186: *>          IU is INTEGER
  187: *>          If RANGE='I', the index of the
  188: *>          largest eigenvalue to be returned.
  189: *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
  190: *>          Not referenced if RANGE = 'A' or 'V'.
  191: *> \endverbatim
  192: *>
  193: *> \param[in] ABSTOL
  194: *> \verbatim
  195: *>          ABSTOL is DOUBLE PRECISION
  196: *>          The absolute error tolerance for the eigenvalues.
  197: *>          An approximate eigenvalue is accepted as converged
  198: *>          when it is determined to lie in an interval [a,b]
  199: *>          of width less than or equal to
  200: *>
  201: *>                  ABSTOL + EPS *   max( |a|,|b| ) ,
  202: *>
  203: *>          where EPS is the machine precision.  If ABSTOL is less than
  204: *>          or equal to zero, then  EPS*|T|  will be used in its place,
  205: *>          where |T| is the 1-norm of the tridiagonal matrix obtained
  206: *>          by reducing A to tridiagonal form.
  207: *>
  208: *>          See "Computing Small Singular Values of Bidiagonal Matrices
  209: *>          with Guaranteed High Relative Accuracy," by Demmel and
  210: *>          Kahan, LAPACK Working Note #3.
  211: *>
  212: *>          If high relative accuracy is important, set ABSTOL to
  213: *>          DLAMCH( 'Safe minimum' ).  Doing so will guarantee that
  214: *>          eigenvalues are computed to high relative accuracy when
  215: *>          possible in future releases.  The current code does not
  216: *>          make any guarantees about high relative accuracy, but
  217: *>          future releases will. See J. Barlow and J. Demmel,
  218: *>          "Computing Accurate Eigensystems of Scaled Diagonally
  219: *>          Dominant Matrices", LAPACK Working Note #7, for a discussion
  220: *>          of which matrices define their eigenvalues to high relative
  221: *>          accuracy.
  222: *> \endverbatim
  223: *>
  224: *> \param[out] M
  225: *> \verbatim
  226: *>          M is INTEGER
  227: *>          The total number of eigenvalues found.  0 <= M <= N.
  228: *>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
  229: *> \endverbatim
  230: *>
  231: *> \param[out] W
  232: *> \verbatim
  233: *>          W is DOUBLE PRECISION array, dimension (N)
  234: *>          The first M elements contain the selected eigenvalues in
  235: *>          ascending order.
  236: *> \endverbatim
  237: *>
  238: *> \param[out] Z
  239: *> \verbatim
  240: *>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
  241: *>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
  242: *>          contain the orthonormal eigenvectors of the matrix A
  243: *>          corresponding to the selected eigenvalues, with the i-th
  244: *>          column of Z holding the eigenvector associated with W(i).
  245: *>          If JOBZ = 'N', then Z is not referenced.
  246: *>          Note: the user must ensure that at least max(1,M) columns are
  247: *>          supplied in the array Z; if RANGE = 'V', the exact value of M
  248: *>          is not known in advance and an upper bound must be used.
  249: *>          Supplying N columns is always safe.
  250: *> \endverbatim
  251: *>
  252: *> \param[in] LDZ
  253: *> \verbatim
  254: *>          LDZ is INTEGER
  255: *>          The leading dimension of the array Z.  LDZ >= 1, and if
  256: *>          JOBZ = 'V', LDZ >= max(1,N).
  257: *> \endverbatim
  258: *>
  259: *> \param[out] ISUPPZ
  260: *> \verbatim
  261: *>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
  262: *>          The support of the eigenvectors in Z, i.e., the indices
  263: *>          indicating the nonzero elements in Z. The i-th eigenvector
  264: *>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
  265: *>          ISUPPZ( 2*i ). This is an output of DSTEMR (tridiagonal
  266: *>          matrix). The support of the eigenvectors of A is typically 
  267: *>          1:N because of the orthogonal transformations applied by DORMTR.
  268: *>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
  269: *> \endverbatim
  270: *>
  271: *> \param[out] WORK
  272: *> \verbatim
  273: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  274: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  275: *> \endverbatim
  276: *>
  277: *> \param[in] LWORK
  278: *> \verbatim
  279: *>          LWORK is INTEGER
  280: *>          The dimension of the array WORK.  
  281: *>          If JOBZ = 'N' and N > 1, LWORK must be queried.
  282: *>                                   LWORK = MAX(1, 26*N, dimension) where
  283: *>                                   dimension = max(stage1,stage2) + (KD+1)*N + 5*N
  284: *>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
  285: *>                                               + max(2*KD*KD, KD*NTHREADS) 
  286: *>                                               + (KD+1)*N + 5*N
  287: *>                                   where KD is the blocking size of the reduction,
  288: *>                                   FACTOPTNB is the blocking used by the QR or LQ
  289: *>                                   algorithm, usually FACTOPTNB=128 is a good choice
  290: *>                                   NTHREADS is the number of threads used when
  291: *>                                   openMP compilation is enabled, otherwise =1.
  292: *>          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available
  293: *>
  294: *>          If LWORK = -1, then a workspace query is assumed; the routine
  295: *>          only calculates the optimal size of the WORK array, returns
  296: *>          this value as the first entry of the WORK array, and no error
  297: *>          message related to LWORK is issued by XERBLA.
  298: *> \endverbatim
  299: *>
  300: *> \param[out] IWORK
  301: *> \verbatim
  302: *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
  303: *>          On exit, if INFO = 0, IWORK(1) returns the optimal LWORK.
  304: *> \endverbatim
  305: *>
  306: *> \param[in] LIWORK
  307: *> \verbatim
  308: *>          LIWORK is INTEGER
  309: *>          The dimension of the array IWORK.  LIWORK >= max(1,10*N).
  310: *>
  311: *>          If LIWORK = -1, then a workspace query is assumed; the
  312: *>          routine only calculates the optimal size of the IWORK array,
  313: *>          returns this value as the first entry of the IWORK array, and
  314: *>          no error message related to LIWORK is issued by XERBLA.
  315: *> \endverbatim
  316: *>
  317: *> \param[out] INFO
  318: *> \verbatim
  319: *>          INFO is INTEGER
  320: *>          = 0:  successful exit
  321: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  322: *>          > 0:  Internal error
  323: *> \endverbatim
  324: *
  325: *  Authors:
  326: *  ========
  327: *
  328: *> \author Univ. of Tennessee
  329: *> \author Univ. of California Berkeley
  330: *> \author Univ. of Colorado Denver
  331: *> \author NAG Ltd.
  332: *
  333: *> \date June 2016
  334: *
  335: *> \ingroup doubleSYeigen
  336: *
  337: *> \par Contributors:
  338: *  ==================
  339: *>
  340: *>     Inderjit Dhillon, IBM Almaden, USA \n
  341: *>     Osni Marques, LBNL/NERSC, USA \n
  342: *>     Ken Stanley, Computer Science Division, University of
  343: *>       California at Berkeley, USA \n
  344: *>     Jason Riedy, Computer Science Division, University of
  345: *>       California at Berkeley, USA \n
  346: *>
  347: *> \par Further Details:
  348: *  =====================
  349: *>
  350: *> \verbatim
  351: *>
  352: *>  All details about the 2stage techniques are available in:
  353: *>
  354: *>  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
  355: *>  Parallel reduction to condensed forms for symmetric eigenvalue problems
  356: *>  using aggregated fine-grained and memory-aware kernels. In Proceedings
  357: *>  of 2011 International Conference for High Performance Computing,
  358: *>  Networking, Storage and Analysis (SC '11), New York, NY, USA,
  359: *>  Article 8 , 11 pages.
  360: *>  http://doi.acm.org/10.1145/2063384.2063394
  361: *>
  362: *>  A. Haidar, J. Kurzak, P. Luszczek, 2013.
  363: *>  An improved parallel singular value algorithm and its implementation 
  364: *>  for multicore hardware, In Proceedings of 2013 International Conference
  365: *>  for High Performance Computing, Networking, Storage and Analysis (SC '13).
  366: *>  Denver, Colorado, USA, 2013.
  367: *>  Article 90, 12 pages.
  368: *>  http://doi.acm.org/10.1145/2503210.2503292
  369: *>
  370: *>  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
  371: *>  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
  372: *>  calculations based on fine-grained memory aware tasks.
  373: *>  International Journal of High Performance Computing Applications.
  374: *>  Volume 28 Issue 2, Pages 196-209, May 2014.
  375: *>  http://hpc.sagepub.com/content/28/2/196 
  376: *>
  377: *> \endverbatim
  378: *
  379: *  =====================================================================
  380:       SUBROUTINE DSYEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
  381:      $                   IL, IU, ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK,
  382:      $                   LWORK, IWORK, LIWORK, INFO )
  383: *
  384:       IMPLICIT NONE
  385: *
  386: *  -- LAPACK driver routine (version 3.7.0) --
  387: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  388: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  389: *     June 2016
  390: *
  391: *     .. Scalar Arguments ..
  392:       CHARACTER          JOBZ, RANGE, UPLO
  393:       INTEGER            IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
  394:       DOUBLE PRECISION   ABSTOL, VL, VU
  395: *     ..
  396: *     .. Array Arguments ..
  397:       INTEGER            ISUPPZ( * ), IWORK( * )
  398:       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
  399: *     ..
  400: *
  401: * =====================================================================
  402: *
  403: *     .. Parameters ..
  404:       DOUBLE PRECISION   ZERO, ONE, TWO
  405:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
  406: *     ..
  407: *     .. Local Scalars ..
  408:       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, VALEIG, WANTZ,
  409:      $                   TRYRAC
  410:       CHARACTER          ORDER
  411:       INTEGER            I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
  412:      $                   INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
  413:      $                   INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
  414:      $                   LLWORK, LLWRKN, LWMIN, NSPLIT,
  415:      $                   LHTRD, LWTRD, KD, IB, INDHOUS
  416:       DOUBLE PRECISION   ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
  417:      $                   SIGMA, SMLNUM, TMP1, VLL, VUU
  418: *     ..
  419: *     .. External Functions ..
  420:       LOGICAL            LSAME
  421:       INTEGER            ILAENV
  422:       DOUBLE PRECISION   DLAMCH, DLANSY
  423:       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANSY
  424: *     ..
  425: *     .. External Subroutines ..
  426:       EXTERNAL           DCOPY, DORMTR, DSCAL, DSTEBZ, DSTEMR, DSTEIN,
  427:      $                   DSTERF, DSWAP, DSYTRD_2STAGE, XERBLA
  428: *     ..
  429: *     .. Intrinsic Functions ..
  430:       INTRINSIC          MAX, MIN, SQRT
  431: *     ..
  432: *     .. Executable Statements ..
  433: *
  434: *     Test the input parameters.
  435: *
  436:       IEEEOK = ILAENV( 10, 'DSYEVR', 'N', 1, 2, 3, 4 )
  437: *
  438:       LOWER = LSAME( UPLO, 'L' )
  439:       WANTZ = LSAME( JOBZ, 'V' )
  440:       ALLEIG = LSAME( RANGE, 'A' )
  441:       VALEIG = LSAME( RANGE, 'V' )
  442:       INDEIG = LSAME( RANGE, 'I' )
  443: *
  444:       LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) )
  445: *
  446:       KD     = ILAENV( 17, 'DSYTRD_2STAGE', JOBZ, N, -1, -1, -1 )
  447:       IB     = ILAENV( 18, 'DSYTRD_2STAGE', JOBZ, N, KD, -1, -1 )
  448:       LHTRD  = ILAENV( 19, 'DSYTRD_2STAGE', JOBZ, N, KD, IB, -1 )
  449:       LWTRD  = ILAENV( 20, 'DSYTRD_2STAGE', JOBZ, N, KD, IB, -1 )
  450:       LWMIN  = MAX( 26*N, 5*N + LHTRD + LWTRD )
  451:       LIWMIN = MAX( 1, 10*N )
  452: *
  453:       INFO = 0
  454:       IF( .NOT.( LSAME( JOBZ, 'N' ) ) ) THEN
  455:          INFO = -1
  456:       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
  457:          INFO = -2
  458:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
  459:          INFO = -3
  460:       ELSE IF( N.LT.0 ) THEN
  461:          INFO = -4
  462:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  463:          INFO = -6
  464:       ELSE
  465:          IF( VALEIG ) THEN
  466:             IF( N.GT.0 .AND. VU.LE.VL )
  467:      $         INFO = -8
  468:          ELSE IF( INDEIG ) THEN
  469:             IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
  470:                INFO = -9
  471:             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
  472:                INFO = -10
  473:             END IF
  474:          END IF
  475:       END IF
  476:       IF( INFO.EQ.0 ) THEN
  477:          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
  478:             INFO = -15
  479:          ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  480:             INFO = -18
  481:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
  482:             INFO = -20
  483:          END IF
  484:       END IF
  485: *
  486:       IF( INFO.EQ.0 ) THEN
  487: *         NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
  488: *         NB = MAX( NB, ILAENV( 1, 'DORMTR', UPLO, N, -1, -1, -1 ) )
  489: *         LWKOPT = MAX( ( NB+1 )*N, LWMIN )
  490:          WORK( 1 ) = LWMIN
  491:          IWORK( 1 ) = LIWMIN
  492:       END IF
  493: *
  494:       IF( INFO.NE.0 ) THEN
  495:          CALL XERBLA( 'DSYEVR_2STAGE', -INFO )
  496:          RETURN
  497:       ELSE IF( LQUERY ) THEN
  498:          RETURN
  499:       END IF
  500: *
  501: *     Quick return if possible
  502: *
  503:       M = 0
  504:       IF( N.EQ.0 ) THEN
  505:          WORK( 1 ) = 1
  506:          RETURN
  507:       END IF
  508: *
  509:       IF( N.EQ.1 ) THEN
  510:          WORK( 1 ) = 7
  511:          IF( ALLEIG .OR. INDEIG ) THEN
  512:             M = 1
  513:             W( 1 ) = A( 1, 1 )
  514:          ELSE
  515:             IF( VL.LT.A( 1, 1 ) .AND. VU.GE.A( 1, 1 ) ) THEN
  516:                M = 1
  517:                W( 1 ) = A( 1, 1 )
  518:             END IF
  519:          END IF
  520:          IF( WANTZ ) THEN
  521:             Z( 1, 1 ) = ONE
  522:             ISUPPZ( 1 ) = 1
  523:             ISUPPZ( 2 ) = 1
  524:          END IF
  525:          RETURN
  526:       END IF
  527: *
  528: *     Get machine constants.
  529: *
  530:       SAFMIN = DLAMCH( 'Safe minimum' )
  531:       EPS    = DLAMCH( 'Precision' )
  532:       SMLNUM = SAFMIN / EPS
  533:       BIGNUM = ONE / SMLNUM
  534:       RMIN   = SQRT( SMLNUM )
  535:       RMAX   = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
  536: *
  537: *     Scale matrix to allowable range, if necessary.
  538: *
  539:       ISCALE = 0
  540:       ABSTLL = ABSTOL
  541:       IF (VALEIG) THEN
  542:          VLL = VL
  543:          VUU = VU
  544:       END IF
  545:       ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
  546:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
  547:          ISCALE = 1
  548:          SIGMA = RMIN / ANRM
  549:       ELSE IF( ANRM.GT.RMAX ) THEN
  550:          ISCALE = 1
  551:          SIGMA = RMAX / ANRM
  552:       END IF
  553:       IF( ISCALE.EQ.1 ) THEN
  554:          IF( LOWER ) THEN
  555:             DO 10 J = 1, N
  556:                CALL DSCAL( N-J+1, SIGMA, A( J, J ), 1 )
  557:    10       CONTINUE
  558:          ELSE
  559:             DO 20 J = 1, N
  560:                CALL DSCAL( J, SIGMA, A( 1, J ), 1 )
  561:    20       CONTINUE
  562:          END IF
  563:          IF( ABSTOL.GT.0 )
  564:      $      ABSTLL = ABSTOL*SIGMA
  565:          IF( VALEIG ) THEN
  566:             VLL = VL*SIGMA
  567:             VUU = VU*SIGMA
  568:          END IF
  569:       END IF
  570: 
  571: *     Initialize indices into workspaces.  Note: The IWORK indices are
  572: *     used only if DSTERF or DSTEMR fail.
  573: 
  574: *     WORK(INDTAU:INDTAU+N-1) stores the scalar factors of the
  575: *     elementary reflectors used in DSYTRD.
  576:       INDTAU = 1
  577: *     WORK(INDD:INDD+N-1) stores the tridiagonal's diagonal entries.
  578:       INDD = INDTAU + N
  579: *     WORK(INDE:INDE+N-1) stores the off-diagonal entries of the
  580: *     tridiagonal matrix from DSYTRD.
  581:       INDE = INDD + N
  582: *     WORK(INDDD:INDDD+N-1) is a copy of the diagonal entries over
  583: *     -written by DSTEMR (the DSTERF path copies the diagonal to W).
  584:       INDDD = INDE + N
  585: *     WORK(INDEE:INDEE+N-1) is a copy of the off-diagonal entries over
  586: *     -written while computing the eigenvalues in DSTERF and DSTEMR.
  587:       INDEE = INDDD + N
  588: *     INDHOUS is the starting offset Householder storage of stage 2
  589:       INDHOUS = INDEE + N
  590: *     INDWK is the starting offset of the left-over workspace, and
  591: *     LLWORK is the remaining workspace size.
  592:       INDWK  = INDHOUS + LHTRD
  593:       LLWORK = LWORK - INDWK + 1
  594: 
  595: 
  596: *     IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
  597: *     stores the block indices of each of the M<=N eigenvalues.
  598:       INDIBL = 1
  599: *     IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
  600: *     stores the starting and finishing indices of each block.
  601:       INDISP = INDIBL + N
  602: *     IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
  603: *     that corresponding to eigenvectors that fail to converge in
  604: *     DSTEIN.  This information is discarded; if any fail, the driver
  605: *     returns INFO > 0.
  606:       INDIFL = INDISP + N
  607: *     INDIWO is the offset of the remaining integer workspace.
  608:       INDIWO = INDIFL + N
  609: 
  610: *
  611: *     Call DSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
  612: *
  613: *
  614:       CALL DSYTRD_2STAGE( JOBZ, UPLO, N, A, LDA, WORK( INDD ), 
  615:      $                    WORK( INDE ), WORK( INDTAU ), WORK( INDHOUS ),
  616:      $                    LHTRD, WORK( INDWK ), LLWORK, IINFO )
  617: *
  618: *     If all eigenvalues are desired
  619: *     then call DSTERF or DSTEMR and DORMTR.
  620: *
  621:       IF( ( ALLEIG .OR. ( INDEIG .AND. IL.EQ.1 .AND. IU.EQ.N ) ) .AND.
  622:      $    IEEEOK.EQ.1 ) THEN
  623:          IF( .NOT.WANTZ ) THEN
  624:             CALL DCOPY( N, WORK( INDD ), 1, W, 1 )
  625:             CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
  626:             CALL DSTERF( N, W, WORK( INDEE ), INFO )
  627:          ELSE
  628:             CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
  629:             CALL DCOPY( N, WORK( INDD ), 1, WORK( INDDD ), 1 )
  630: *
  631:             IF (ABSTOL .LE. TWO*N*EPS) THEN
  632:                TRYRAC = .TRUE.
  633:             ELSE
  634:                TRYRAC = .FALSE.
  635:             END IF
  636:             CALL DSTEMR( JOBZ, 'A', N, WORK( INDDD ), WORK( INDEE ),
  637:      $                   VL, VU, IL, IU, M, W, Z, LDZ, N, ISUPPZ,
  638:      $                   TRYRAC, WORK( INDWK ), LWORK, IWORK, LIWORK,
  639:      $                   INFO )
  640: *
  641: *
  642: *
  643: *        Apply orthogonal matrix used in reduction to tridiagonal
  644: *        form to eigenvectors returned by DSTEMR.
  645: *
  646:             IF( WANTZ .AND. INFO.EQ.0 ) THEN
  647:                INDWKN = INDE
  648:                LLWRKN = LWORK - INDWKN + 1
  649:                CALL DORMTR( 'L', UPLO, 'N', N, M, A, LDA,
  650:      $                      WORK( INDTAU ), Z, LDZ, WORK( INDWKN ),
  651:      $                      LLWRKN, IINFO )
  652:             END IF
  653:          END IF
  654: *
  655: *
  656:          IF( INFO.EQ.0 ) THEN
  657: *           Everything worked.  Skip DSTEBZ/DSTEIN.  IWORK(:) are
  658: *           undefined.
  659:             M = N
  660:             GO TO 30
  661:          END IF
  662:          INFO = 0
  663:       END IF
  664: *
  665: *     Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN.
  666: *     Also call DSTEBZ and DSTEIN if DSTEMR fails.
  667: *
  668:       IF( WANTZ ) THEN
  669:          ORDER = 'B'
  670:       ELSE
  671:          ORDER = 'E'
  672:       END IF
  673: 
  674:       CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
  675:      $             WORK( INDD ), WORK( INDE ), M, NSPLIT, W,
  676:      $             IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWK ),
  677:      $             IWORK( INDIWO ), INFO )
  678: *
  679:       IF( WANTZ ) THEN
  680:          CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W,
  681:      $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
  682:      $                WORK( INDWK ), IWORK( INDIWO ), IWORK( INDIFL ),
  683:      $                INFO )
  684: *
  685: *        Apply orthogonal matrix used in reduction to tridiagonal
  686: *        form to eigenvectors returned by DSTEIN.
  687: *
  688:          INDWKN = INDE
  689:          LLWRKN = LWORK - INDWKN + 1
  690:          CALL DORMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z,
  691:      $                LDZ, WORK( INDWKN ), LLWRKN, IINFO )
  692:       END IF
  693: *
  694: *     If matrix was scaled, then rescale eigenvalues appropriately.
  695: *
  696: *  Jump here if DSTEMR/DSTEIN succeeded.
  697:    30 CONTINUE
  698:       IF( ISCALE.EQ.1 ) THEN
  699:          IF( INFO.EQ.0 ) THEN
  700:             IMAX = M
  701:          ELSE
  702:             IMAX = INFO - 1
  703:          END IF
  704:          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
  705:       END IF
  706: *
  707: *     If eigenvalues are not in order, then sort them, along with
  708: *     eigenvectors.  Note: We do not sort the IFAIL portion of IWORK.
  709: *     It may not be initialized (if DSTEMR/DSTEIN succeeded), and we do
  710: *     not return this detailed information to the user.
  711: *
  712:       IF( WANTZ ) THEN
  713:          DO 50 J = 1, M - 1
  714:             I = 0
  715:             TMP1 = W( J )
  716:             DO 40 JJ = J + 1, M
  717:                IF( W( JJ ).LT.TMP1 ) THEN
  718:                   I = JJ
  719:                   TMP1 = W( JJ )
  720:                END IF
  721:    40       CONTINUE
  722: *
  723:             IF( I.NE.0 ) THEN
  724:                W( I ) = W( J )
  725:                W( J ) = TMP1
  726:                CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
  727:             END IF
  728:    50    CONTINUE
  729:       END IF
  730: *
  731: *     Set WORK(1) to optimal workspace size.
  732: *
  733:       WORK( 1 ) = LWMIN
  734:       IWORK( 1 ) = LIWMIN
  735: *
  736:       RETURN
  737: *
  738: *     End of DSYEVR_2STAGE
  739: *
  740:       END

CVSweb interface <joel.bertrand@systella.fr>