Annotation of rpl/lapack/lapack/zhbgvd.f, revision 1.9

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

CVSweb interface <joel.bertrand@systella.fr>