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

1.1     ! bertrand    1: *> \brief <b> DSYEV_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 DSYEV_2STAGE + dependencies
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyevd_2stage.f">
        !            13: *> [TGZ]</a>
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyevd_2stage.f">
        !            15: *> [ZIP]</a>
        !            16: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyevd_2stage.f">
        !            17: *> [TXT]</a>
        !            18: *> \endhtmlonly
        !            19: *
        !            20: *  Definition:
        !            21: *  ===========
        !            22: *
        !            23: *       SUBROUTINE DSYEV_2STAGE( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, 
        !            24: *                                INFO )
        !            25: *
        !            26: *       IMPLICIT NONE
        !            27: *
        !            28: *       .. Scalar Arguments ..
        !            29: *       CHARACTER          JOBZ, UPLO
        !            30: *       INTEGER            INFO, LDA, LWORK, N
        !            31: *       ..
        !            32: *       .. Array Arguments ..
        !            33: *       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
        !            34: *       ..
        !            35: *
        !            36: *
        !            37: *> \par Purpose:
        !            38: *  =============
        !            39: *>
        !            40: *> \verbatim
        !            41: *>
        !            42: *> DSYEV_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
        !            43: *> real symmetric matrix A using the 2stage technique for
        !            44: *> the reduction to tridiagonal.
        !            45: *> \endverbatim
        !            46: *
        !            47: *  Arguments:
        !            48: *  ==========
        !            49: *
        !            50: *> \param[in] JOBZ
        !            51: *> \verbatim
        !            52: *>          JOBZ is CHARACTER*1
        !            53: *>          = 'N':  Compute eigenvalues only;
        !            54: *>          = 'V':  Compute eigenvalues and eigenvectors.
        !            55: *>                  Not available in this release.
        !            56: *> \endverbatim
        !            57: *>
        !            58: *> \param[in] UPLO
        !            59: *> \verbatim
        !            60: *>          UPLO is CHARACTER*1
        !            61: *>          = 'U':  Upper triangle of A is stored;
        !            62: *>          = 'L':  Lower triangle of A is stored.
        !            63: *> \endverbatim
        !            64: *>
        !            65: *> \param[in] N
        !            66: *> \verbatim
        !            67: *>          N is INTEGER
        !            68: *>          The order of the matrix A.  N >= 0.
        !            69: *> \endverbatim
        !            70: *>
        !            71: *> \param[in,out] A
        !            72: *> \verbatim
        !            73: *>          A is DOUBLE PRECISION array, dimension (LDA, N)
        !            74: *>          On entry, the symmetric matrix A.  If UPLO = 'U', the
        !            75: *>          leading N-by-N upper triangular part of A contains the
        !            76: *>          upper triangular part of the matrix A.  If UPLO = 'L',
        !            77: *>          the leading N-by-N lower triangular part of A contains
        !            78: *>          the lower triangular part of the matrix A.
        !            79: *>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
        !            80: *>          orthonormal eigenvectors of the matrix A.
        !            81: *>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
        !            82: *>          or the upper triangle (if UPLO='U') of A, including the
        !            83: *>          diagonal, is destroyed.
        !            84: *> \endverbatim
        !            85: *>
        !            86: *> \param[in] LDA
        !            87: *> \verbatim
        !            88: *>          LDA is INTEGER
        !            89: *>          The leading dimension of the array A.  LDA >= max(1,N).
        !            90: *> \endverbatim
        !            91: *>
        !            92: *> \param[out] W
        !            93: *> \verbatim
        !            94: *>          W is DOUBLE PRECISION array, dimension (N)
        !            95: *>          If INFO = 0, the eigenvalues in ascending order.
        !            96: *> \endverbatim
        !            97: *>
        !            98: *> \param[out] WORK
        !            99: *> \verbatim
        !           100: *>          WORK is DOUBLE PRECISION array, dimension LWORK
        !           101: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
        !           102: *> \endverbatim
        !           103: *>
        !           104: *> \param[in] LWORK
        !           105: *> \verbatim
        !           106: *>          LWORK is INTEGER
        !           107: *>          The length of the array WORK. LWORK >= 1, when N <= 1;
        !           108: *>          otherwise  
        !           109: *>          If JOBZ = 'N' and N > 1, LWORK must be queried.
        !           110: *>                                   LWORK = MAX(1, dimension) where
        !           111: *>                                   dimension = max(stage1,stage2) + (KD+1)*N + 2*N
        !           112: *>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
        !           113: *>                                               + max(2*KD*KD, KD*NTHREADS) 
        !           114: *>                                               + (KD+1)*N + 2*N
        !           115: *>                                   where KD is the blocking size of the reduction,
        !           116: *>                                   FACTOPTNB is the blocking used by the QR or LQ
        !           117: *>                                   algorithm, usually FACTOPTNB=128 is a good choice
        !           118: *>                                   NTHREADS is the number of threads used when
        !           119: *>                                   openMP compilation is enabled, otherwise =1.
        !           120: *>          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available
        !           121: *>
        !           122: *>          If LWORK = -1, then a workspace query is assumed; the routine
        !           123: *>          only calculates the optimal size of the WORK array, returns
        !           124: *>          this value as the first entry of the WORK array, and no error
        !           125: *>          message related to LWORK is issued by XERBLA.
        !           126: *> \endverbatim
        !           127: *>
        !           128: *> \param[out] INFO
        !           129: *> \verbatim
        !           130: *>          INFO is INTEGER
        !           131: *>          = 0:  successful exit
        !           132: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
        !           133: *>          > 0:  if INFO = i, the algorithm failed to converge; i
        !           134: *>                off-diagonal elements of an intermediate tridiagonal
        !           135: *>                form did not converge to zero.
        !           136: *> \endverbatim
        !           137: *
        !           138: *  Authors:
        !           139: *  ========
        !           140: *
        !           141: *> \author Univ. of Tennessee
        !           142: *> \author Univ. of California Berkeley
        !           143: *> \author Univ. of Colorado Denver
        !           144: *> \author NAG Ltd.
        !           145: *
        !           146: *> \date December 2016
        !           147: *
        !           148: *> \ingroup doubleSYeigen
        !           149: *
        !           150: *> \par Further Details:
        !           151: *  =====================
        !           152: *>
        !           153: *> \verbatim
        !           154: *>
        !           155: *>  All details about the 2stage techniques are available in:
        !           156: *>
        !           157: *>  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
        !           158: *>  Parallel reduction to condensed forms for symmetric eigenvalue problems
        !           159: *>  using aggregated fine-grained and memory-aware kernels. In Proceedings
        !           160: *>  of 2011 International Conference for High Performance Computing,
        !           161: *>  Networking, Storage and Analysis (SC '11), New York, NY, USA,
        !           162: *>  Article 8 , 11 pages.
        !           163: *>  http://doi.acm.org/10.1145/2063384.2063394
        !           164: *>
        !           165: *>  A. Haidar, J. Kurzak, P. Luszczek, 2013.
        !           166: *>  An improved parallel singular value algorithm and its implementation 
        !           167: *>  for multicore hardware, In Proceedings of 2013 International Conference
        !           168: *>  for High Performance Computing, Networking, Storage and Analysis (SC '13).
        !           169: *>  Denver, Colorado, USA, 2013.
        !           170: *>  Article 90, 12 pages.
        !           171: *>  http://doi.acm.org/10.1145/2503210.2503292
        !           172: *>
        !           173: *>  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
        !           174: *>  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
        !           175: *>  calculations based on fine-grained memory aware tasks.
        !           176: *>  International Journal of High Performance Computing Applications.
        !           177: *>  Volume 28 Issue 2, Pages 196-209, May 2014.
        !           178: *>  http://hpc.sagepub.com/content/28/2/196 
        !           179: *>
        !           180: *> \endverbatim
        !           181: *
        !           182: *  =====================================================================
        !           183:       SUBROUTINE DSYEV_2STAGE( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, 
        !           184:      $                         INFO )
        !           185: *
        !           186:       IMPLICIT NONE
        !           187: *
        !           188: *  -- LAPACK driver routine (version 3.7.0) --
        !           189: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
        !           190: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
        !           191: *     December 2016
        !           192: *
        !           193: *     .. Scalar Arguments ..
        !           194:       CHARACTER          JOBZ, UPLO
        !           195:       INTEGER            INFO, LDA, LWORK, N
        !           196: *     ..
        !           197: *     .. Array Arguments ..
        !           198:       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
        !           199: *     ..
        !           200: *
        !           201: *  =====================================================================
        !           202: *
        !           203: *     .. Parameters ..
        !           204:       DOUBLE PRECISION   ZERO, ONE
        !           205:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
        !           206: *     ..
        !           207: *     .. Local Scalars ..
        !           208:       LOGICAL            LOWER, LQUERY, WANTZ
        !           209:       INTEGER            IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
        !           210:      $                   LLWORK, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
        !           211:       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
        !           212:      $                   SMLNUM
        !           213: *     ..
        !           214: *     .. External Functions ..
        !           215:       LOGICAL            LSAME
        !           216:       INTEGER            ILAENV
        !           217:       DOUBLE PRECISION   DLAMCH, DLANSY
        !           218:       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANSY
        !           219: *     ..
        !           220: *     .. External Subroutines ..
        !           221:       EXTERNAL           DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF,
        !           222:      $                   XERBLA, DSYTRD_2STAGE
        !           223: *     ..
        !           224: *     .. Intrinsic Functions ..
        !           225:       INTRINSIC          MAX, SQRT
        !           226: *     ..
        !           227: *     .. Executable Statements ..
        !           228: *
        !           229: *     Test the input parameters.
        !           230: *
        !           231:       WANTZ = LSAME( JOBZ, 'V' )
        !           232:       LOWER = LSAME( UPLO, 'L' )
        !           233:       LQUERY = ( LWORK.EQ.-1 )
        !           234: *
        !           235:       INFO = 0
        !           236:       IF( .NOT.( LSAME( JOBZ, 'N' ) ) ) THEN
        !           237:          INFO = -1
        !           238:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
        !           239:          INFO = -2
        !           240:       ELSE IF( N.LT.0 ) THEN
        !           241:          INFO = -3
        !           242:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
        !           243:          INFO = -5
        !           244:       END IF
        !           245: *
        !           246:       IF( INFO.EQ.0 ) THEN
        !           247:          KD    = ILAENV( 17, 'DSYTRD_2STAGE', JOBZ, N, -1, -1, -1 )
        !           248:          IB    = ILAENV( 18, 'DSYTRD_2STAGE', JOBZ, N, KD, -1, -1 )
        !           249:          LHTRD = ILAENV( 19, 'DSYTRD_2STAGE', JOBZ, N, KD, IB, -1 )
        !           250:          LWTRD = ILAENV( 20, 'DSYTRD_2STAGE', JOBZ, N, KD, IB, -1 )
        !           251:          LWMIN = 2*N + LHTRD + LWTRD
        !           252:          WORK( 1 )  = LWMIN
        !           253: *
        !           254:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY )
        !           255:      $      INFO = -8
        !           256:       END IF
        !           257: *
        !           258:       IF( INFO.NE.0 ) THEN
        !           259:          CALL XERBLA( 'DSYEV_2STAGE ', -INFO )
        !           260:          RETURN
        !           261:       ELSE IF( LQUERY ) THEN
        !           262:          RETURN
        !           263:       END IF
        !           264: *
        !           265: *     Quick return if possible
        !           266: *
        !           267:       IF( N.EQ.0 ) THEN
        !           268:          RETURN
        !           269:       END IF
        !           270: *
        !           271:       IF( N.EQ.1 ) THEN
        !           272:          W( 1 ) = A( 1, 1 )
        !           273:          WORK( 1 ) = 2
        !           274:          IF( WANTZ )
        !           275:      $      A( 1, 1 ) = ONE
        !           276:          RETURN
        !           277:       END IF
        !           278: *
        !           279: *     Get machine constants.
        !           280: *
        !           281:       SAFMIN = DLAMCH( 'Safe minimum' )
        !           282:       EPS    = DLAMCH( 'Precision' )
        !           283:       SMLNUM = SAFMIN / EPS
        !           284:       BIGNUM = ONE / SMLNUM
        !           285:       RMIN   = SQRT( SMLNUM )
        !           286:       RMAX   = SQRT( BIGNUM )
        !           287: *
        !           288: *     Scale matrix to allowable range, if necessary.
        !           289: *
        !           290:       ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
        !           291:       ISCALE = 0
        !           292:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
        !           293:          ISCALE = 1
        !           294:          SIGMA = RMIN / ANRM
        !           295:       ELSE IF( ANRM.GT.RMAX ) THEN
        !           296:          ISCALE = 1
        !           297:          SIGMA = RMAX / ANRM
        !           298:       END IF
        !           299:       IF( ISCALE.EQ.1 )
        !           300:      $   CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
        !           301: *
        !           302: *     Call DSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
        !           303: *
        !           304:       INDE    = 1
        !           305:       INDTAU  = INDE + N
        !           306:       INDHOUS = INDTAU + N
        !           307:       INDWRK  = INDHOUS + LHTRD
        !           308:       LLWORK  = LWORK - INDWRK + 1
        !           309: *
        !           310:       CALL DSYTRD_2STAGE( JOBZ, UPLO, N, A, LDA, W, WORK( INDE ),
        !           311:      $                    WORK( INDTAU ), WORK( INDHOUS ), LHTRD, 
        !           312:      $                    WORK( INDWRK ), LLWORK, IINFO )
        !           313: *
        !           314: *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
        !           315: *     DORGTR to generate the orthogonal matrix, then call DSTEQR.
        !           316: *
        !           317:       IF( .NOT.WANTZ ) THEN
        !           318:          CALL DSTERF( N, W, WORK( INDE ), INFO )
        !           319:       ELSE
        !           320: *        Not available in this release, and agrument checking should not
        !           321: *        let it getting here
        !           322:          RETURN
        !           323:          CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
        !           324:      $                LLWORK, IINFO )
        !           325:          CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ),
        !           326:      $                INFO )
        !           327:       END IF
        !           328: *
        !           329: *     If matrix was scaled, then rescale eigenvalues appropriately.
        !           330: *
        !           331:       IF( ISCALE.EQ.1 ) THEN
        !           332:          IF( INFO.EQ.0 ) THEN
        !           333:             IMAX = N
        !           334:          ELSE
        !           335:             IMAX = INFO - 1
        !           336:          END IF
        !           337:          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
        !           338:       END IF
        !           339: *
        !           340: *     Set WORK(1) to optimal workspace size.
        !           341: *
        !           342:       WORK( 1 ) = LWMIN
        !           343: *
        !           344:       RETURN
        !           345: *
        !           346: *     End of DSYEV_2STAGE
        !           347: *
        !           348:       END

CVSweb interface <joel.bertrand@systella.fr>