Annotation of rpl/lapack/lapack/zhbevd.f, revision 1.8

1.8     ! bertrand    1: *> \brief <b> ZHBEVD 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 ZHBEVD + dependencies 
        !            10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbevd.f"> 
        !            11: *> [TGZ]</a> 
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbevd.f"> 
        !            13: *> [ZIP]</a> 
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbevd.f"> 
        !            15: *> [TXT]</a>
        !            16: *> \endhtmlonly 
        !            17: *
        !            18: *  Definition:
        !            19: *  ===========
        !            20: *
        !            21: *       SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
        !            22: *                          LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
        !            23: * 
        !            24: *       .. Scalar Arguments ..
        !            25: *       CHARACTER          JOBZ, UPLO
        !            26: *       INTEGER            INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
        !            27: *       ..
        !            28: *       .. Array Arguments ..
        !            29: *       INTEGER            IWORK( * )
        !            30: *       DOUBLE PRECISION   RWORK( * ), W( * )
        !            31: *       COMPLEX*16         AB( LDAB, * ), WORK( * ), Z( LDZ, * )
        !            32: *       ..
        !            33: *  
        !            34: *
        !            35: *> \par Purpose:
        !            36: *  =============
        !            37: *>
        !            38: *> \verbatim
        !            39: *>
        !            40: *> ZHBEVD computes all the eigenvalues and, optionally, eigenvectors of
        !            41: *> a complex Hermitian band matrix A.  If eigenvectors are desired, it
        !            42: *> uses a divide and conquer algorithm.
        !            43: *>
        !            44: *> The divide and conquer algorithm makes very mild assumptions about
        !            45: *> floating point arithmetic. It will work on machines with a guard
        !            46: *> digit in add/subtract, or on those binary machines without guard
        !            47: *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
        !            48: *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
        !            49: *> without guard digits, but we know of none.
        !            50: *> \endverbatim
        !            51: *
        !            52: *  Arguments:
        !            53: *  ==========
        !            54: *
        !            55: *> \param[in] JOBZ
        !            56: *> \verbatim
        !            57: *>          JOBZ is CHARACTER*1
        !            58: *>          = 'N':  Compute eigenvalues only;
        !            59: *>          = 'V':  Compute eigenvalues and eigenvectors.
        !            60: *> \endverbatim
        !            61: *>
        !            62: *> \param[in] UPLO
        !            63: *> \verbatim
        !            64: *>          UPLO is CHARACTER*1
        !            65: *>          = 'U':  Upper triangle of A is stored;
        !            66: *>          = 'L':  Lower triangle of A is stored.
        !            67: *> \endverbatim
        !            68: *>
        !            69: *> \param[in] N
        !            70: *> \verbatim
        !            71: *>          N is INTEGER
        !            72: *>          The order of the matrix A.  N >= 0.
        !            73: *> \endverbatim
        !            74: *>
        !            75: *> \param[in] KD
        !            76: *> \verbatim
        !            77: *>          KD is INTEGER
        !            78: *>          The number of superdiagonals of the matrix A if UPLO = 'U',
        !            79: *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
        !            80: *> \endverbatim
        !            81: *>
        !            82: *> \param[in,out] AB
        !            83: *> \verbatim
        !            84: *>          AB is COMPLEX*16 array, dimension (LDAB, N)
        !            85: *>          On entry, the upper or lower triangle of the Hermitian band
        !            86: *>          matrix A, stored in the first KD+1 rows of the array.  The
        !            87: *>          j-th column of A is stored in the j-th column of the array AB
        !            88: *>          as follows:
        !            89: *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
        !            90: *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
        !            91: *>
        !            92: *>          On exit, AB is overwritten by values generated during the
        !            93: *>          reduction to tridiagonal form.  If UPLO = 'U', the first
        !            94: *>          superdiagonal and the diagonal of the tridiagonal matrix T
        !            95: *>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
        !            96: *>          the diagonal and first subdiagonal of T are returned in the
        !            97: *>          first two rows of AB.
        !            98: *> \endverbatim
        !            99: *>
        !           100: *> \param[in] LDAB
        !           101: *> \verbatim
        !           102: *>          LDAB is INTEGER
        !           103: *>          The leading dimension of the array AB.  LDAB >= KD + 1.
        !           104: *> \endverbatim
        !           105: *>
        !           106: *> \param[out] W
        !           107: *> \verbatim
        !           108: *>          W is DOUBLE PRECISION array, dimension (N)
        !           109: *>          If INFO = 0, the eigenvalues in ascending order.
        !           110: *> \endverbatim
        !           111: *>
        !           112: *> \param[out] Z
        !           113: *> \verbatim
        !           114: *>          Z is COMPLEX*16 array, dimension (LDZ, N)
        !           115: *>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
        !           116: *>          eigenvectors of the matrix A, with the i-th column of Z
        !           117: *>          holding the eigenvector associated with W(i).
        !           118: *>          If JOBZ = 'N', then Z is not referenced.
        !           119: *> \endverbatim
        !           120: *>
        !           121: *> \param[in] LDZ
        !           122: *> \verbatim
        !           123: *>          LDZ is INTEGER
        !           124: *>          The leading dimension of the array Z.  LDZ >= 1, and if
        !           125: *>          JOBZ = 'V', LDZ >= max(1,N).
        !           126: *> \endverbatim
        !           127: *>
        !           128: *> \param[out] WORK
        !           129: *> \verbatim
        !           130: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
        !           131: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
        !           132: *> \endverbatim
        !           133: *>
        !           134: *> \param[in] LWORK
        !           135: *> \verbatim
        !           136: *>          LWORK is INTEGER
        !           137: *>          The dimension of the array WORK.
        !           138: *>          If N <= 1,               LWORK must be at least 1.
        !           139: *>          If JOBZ = 'N' and N > 1, LWORK must be at least N.
        !           140: *>          If JOBZ = 'V' and N > 1, LWORK must be at least 2*N**2.
        !           141: *>
        !           142: *>          If LWORK = -1, then a workspace query is assumed; the routine
        !           143: *>          only calculates the optimal sizes of the WORK, RWORK and
        !           144: *>          IWORK arrays, returns these values as the first entries of
        !           145: *>          the WORK, RWORK and IWORK arrays, and no error message
        !           146: *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
        !           147: *> \endverbatim
        !           148: *>
        !           149: *> \param[out] RWORK
        !           150: *> \verbatim
        !           151: *>          RWORK is DOUBLE PRECISION array,
        !           152: *>                                         dimension (LRWORK)
        !           153: *>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
        !           154: *> \endverbatim
        !           155: *>
        !           156: *> \param[in] LRWORK
        !           157: *> \verbatim
        !           158: *>          LRWORK is INTEGER
        !           159: *>          The dimension of array RWORK.
        !           160: *>          If N <= 1,               LRWORK must be at least 1.
        !           161: *>          If JOBZ = 'N' and N > 1, LRWORK must be at least N.
        !           162: *>          If JOBZ = 'V' and N > 1, LRWORK must be at least
        !           163: *>                        1 + 5*N + 2*N**2.
        !           164: *>
        !           165: *>          If LRWORK = -1, then a workspace query is assumed; the
        !           166: *>          routine only calculates the optimal sizes of the WORK, RWORK
        !           167: *>          and IWORK arrays, returns these values as the first entries
        !           168: *>          of the WORK, RWORK and IWORK arrays, and no error message
        !           169: *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
        !           170: *> \endverbatim
        !           171: *>
        !           172: *> \param[out] IWORK
        !           173: *> \verbatim
        !           174: *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
        !           175: *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
        !           176: *> \endverbatim
        !           177: *>
        !           178: *> \param[in] LIWORK
        !           179: *> \verbatim
        !           180: *>          LIWORK is INTEGER
        !           181: *>          The dimension of array IWORK.
        !           182: *>          If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
        !           183: *>          If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N .
        !           184: *>
        !           185: *>          If LIWORK = -1, then a workspace query is assumed; the
        !           186: *>          routine only calculates the optimal sizes of the WORK, RWORK
        !           187: *>          and IWORK arrays, returns these values as the first entries
        !           188: *>          of the WORK, RWORK and IWORK arrays, and no error message
        !           189: *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
        !           190: *> \endverbatim
        !           191: *>
        !           192: *> \param[out] INFO
        !           193: *> \verbatim
        !           194: *>          INFO is INTEGER
        !           195: *>          = 0:  successful exit.
        !           196: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
        !           197: *>          > 0:  if INFO = i, the algorithm failed to converge; i
        !           198: *>                off-diagonal elements of an intermediate tridiagonal
        !           199: *>                form did not converge to zero.
        !           200: *> \endverbatim
        !           201: *
        !           202: *  Authors:
        !           203: *  ========
        !           204: *
        !           205: *> \author Univ. of Tennessee 
        !           206: *> \author Univ. of California Berkeley 
        !           207: *> \author Univ. of Colorado Denver 
        !           208: *> \author NAG Ltd. 
        !           209: *
        !           210: *> \date November 2011
        !           211: *
        !           212: *> \ingroup complex16OTHEReigen
        !           213: *
        !           214: *  =====================================================================
1.1       bertrand  215:       SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
                    216:      $                   LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
                    217: *
1.8     ! bertrand  218: *  -- LAPACK driver routine (version 3.4.0) --
1.1       bertrand  219: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    220: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8     ! bertrand  221: *     November 2011
1.1       bertrand  222: *
                    223: *     .. Scalar Arguments ..
                    224:       CHARACTER          JOBZ, UPLO
                    225:       INTEGER            INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
                    226: *     ..
                    227: *     .. Array Arguments ..
                    228:       INTEGER            IWORK( * )
                    229:       DOUBLE PRECISION   RWORK( * ), W( * )
                    230:       COMPLEX*16         AB( LDAB, * ), WORK( * ), Z( LDZ, * )
                    231: *     ..
                    232: *
                    233: *  =====================================================================
                    234: *
                    235: *     .. Parameters ..
                    236:       DOUBLE PRECISION   ZERO, ONE
                    237:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
                    238:       COMPLEX*16         CZERO, CONE
                    239:       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
                    240:      $                   CONE = ( 1.0D0, 0.0D0 ) )
                    241: *     ..
                    242: *     .. Local Scalars ..
                    243:       LOGICAL            LOWER, LQUERY, WANTZ
                    244:       INTEGER            IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE,
                    245:      $                   LIWMIN, LLRWK, LLWK2, LRWMIN, LWMIN
                    246:       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
                    247:      $                   SMLNUM
                    248: *     ..
                    249: *     .. External Functions ..
                    250:       LOGICAL            LSAME
                    251:       DOUBLE PRECISION   DLAMCH, ZLANHB
                    252:       EXTERNAL           LSAME, DLAMCH, ZLANHB
                    253: *     ..
                    254: *     .. External Subroutines ..
                    255:       EXTERNAL           DSCAL, DSTERF, XERBLA, ZGEMM, ZHBTRD, ZLACPY,
                    256:      $                   ZLASCL, ZSTEDC
                    257: *     ..
                    258: *     .. Intrinsic Functions ..
                    259:       INTRINSIC          SQRT
                    260: *     ..
                    261: *     .. Executable Statements ..
                    262: *
                    263: *     Test the input parameters.
                    264: *
                    265:       WANTZ = LSAME( JOBZ, 'V' )
                    266:       LOWER = LSAME( UPLO, 'L' )
                    267:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
                    268: *
                    269:       INFO = 0
                    270:       IF( N.LE.1 ) THEN
                    271:          LWMIN = 1
                    272:          LRWMIN = 1
                    273:          LIWMIN = 1
                    274:       ELSE
                    275:          IF( WANTZ ) THEN
                    276:             LWMIN = 2*N**2
                    277:             LRWMIN = 1 + 5*N + 2*N**2
                    278:             LIWMIN = 3 + 5*N
                    279:          ELSE
                    280:             LWMIN = N
                    281:             LRWMIN = N
                    282:             LIWMIN = 1
                    283:          END IF
                    284:       END IF
                    285:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
                    286:          INFO = -1
                    287:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
                    288:          INFO = -2
                    289:       ELSE IF( N.LT.0 ) THEN
                    290:          INFO = -3
                    291:       ELSE IF( KD.LT.0 ) THEN
                    292:          INFO = -4
                    293:       ELSE IF( LDAB.LT.KD+1 ) THEN
                    294:          INFO = -6
                    295:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
                    296:          INFO = -9
                    297:       END IF
                    298: *
                    299:       IF( INFO.EQ.0 ) THEN
                    300:          WORK( 1 ) = LWMIN
                    301:          RWORK( 1 ) = LRWMIN
                    302:          IWORK( 1 ) = LIWMIN
                    303: *
                    304:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
                    305:             INFO = -11
                    306:          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
                    307:             INFO = -13
                    308:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
                    309:             INFO = -15
                    310:          END IF
                    311:       END IF
                    312: *
                    313:       IF( INFO.NE.0 ) THEN
                    314:          CALL XERBLA( 'ZHBEVD', -INFO )
                    315:          RETURN
                    316:       ELSE IF( LQUERY ) THEN
                    317:          RETURN
                    318:       END IF
                    319: *
                    320: *     Quick return if possible
                    321: *
                    322:       IF( N.EQ.0 )
                    323:      $   RETURN
                    324: *
                    325:       IF( N.EQ.1 ) THEN
                    326:          W( 1 ) = AB( 1, 1 )
                    327:          IF( WANTZ )
                    328:      $      Z( 1, 1 ) = CONE
                    329:          RETURN
                    330:       END IF
                    331: *
                    332: *     Get machine constants.
                    333: *
                    334:       SAFMIN = DLAMCH( 'Safe minimum' )
                    335:       EPS = DLAMCH( 'Precision' )
                    336:       SMLNUM = SAFMIN / EPS
                    337:       BIGNUM = ONE / SMLNUM
                    338:       RMIN = SQRT( SMLNUM )
                    339:       RMAX = SQRT( BIGNUM )
                    340: *
                    341: *     Scale matrix to allowable range, if necessary.
                    342: *
                    343:       ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
                    344:       ISCALE = 0
                    345:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
                    346:          ISCALE = 1
                    347:          SIGMA = RMIN / ANRM
                    348:       ELSE IF( ANRM.GT.RMAX ) THEN
                    349:          ISCALE = 1
                    350:          SIGMA = RMAX / ANRM
                    351:       END IF
                    352:       IF( ISCALE.EQ.1 ) THEN
                    353:          IF( LOWER ) THEN
                    354:             CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
                    355:          ELSE
                    356:             CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
                    357:          END IF
                    358:       END IF
                    359: *
                    360: *     Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
                    361: *
                    362:       INDE = 1
                    363:       INDWRK = INDE + N
                    364:       INDWK2 = 1 + N*N
                    365:       LLWK2 = LWORK - INDWK2 + 1
                    366:       LLRWK = LRWORK - INDWRK + 1
                    367:       CALL ZHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z,
                    368:      $             LDZ, WORK, IINFO )
                    369: *
                    370: *     For eigenvalues only, call DSTERF.  For eigenvectors, call ZSTEDC.
                    371: *
                    372:       IF( .NOT.WANTZ ) THEN
                    373:          CALL DSTERF( N, W, RWORK( INDE ), INFO )
                    374:       ELSE
                    375:          CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
                    376:      $                LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
                    377:      $                INFO )
                    378:          CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
                    379:      $               WORK( INDWK2 ), N )
                    380:          CALL ZLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
                    381:       END IF
                    382: *
                    383: *     If matrix was scaled, then rescale eigenvalues appropriately.
                    384: *
                    385:       IF( ISCALE.EQ.1 ) THEN
                    386:          IF( INFO.EQ.0 ) THEN
                    387:             IMAX = N
                    388:          ELSE
                    389:             IMAX = INFO - 1
                    390:          END IF
                    391:          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
                    392:       END IF
                    393: *
                    394:       WORK( 1 ) = LWMIN
                    395:       RWORK( 1 ) = LRWMIN
                    396:       IWORK( 1 ) = LIWMIN
                    397:       RETURN
                    398: *
                    399: *     End of ZHBEVD
                    400: *
                    401:       END

CVSweb interface <joel.bertrand@systella.fr>