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

1.14      bertrand    1: *> \brief \b ZHBGVD
1.9       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.17      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.9       bertrand    7: *
                      8: *> \htmlonly
1.17      bertrand    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">
1.9       bertrand   15: *> [TXT]</a>
1.17      bertrand   16: *> \endhtmlonly
1.9       bertrand   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 )
1.17      bertrand   24: *
1.9       bertrand   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: *       ..
1.17      bertrand   36: *
1.9       bertrand   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: *
1.17      bertrand  236: *> \author Univ. of Tennessee
                    237: *> \author Univ. of California Berkeley
                    238: *> \author Univ. of Colorado Denver
                    239: *> \author NAG Ltd.
1.9       bertrand  240: *
                    241: *> \ingroup complex16OTHEReigen
                    242: *
                    243: *> \par Contributors:
                    244: *  ==================
                    245: *>
                    246: *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
                    247: *
                    248: *  =====================================================================
1.1       bertrand  249:       SUBROUTINE ZHBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
                    250:      $                   Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK,
                    251:      $                   LIWORK, INFO )
                    252: *
1.20    ! bertrand  253: *  -- LAPACK driver routine --
1.1       bertrand  254: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    255: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    256: *
                    257: *     .. Scalar Arguments ..
                    258:       CHARACTER          JOBZ, UPLO
                    259:       INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LRWORK,
                    260:      $                   LWORK, N
                    261: *     ..
                    262: *     .. Array Arguments ..
                    263:       INTEGER            IWORK( * )
                    264:       DOUBLE PRECISION   RWORK( * ), W( * )
                    265:       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
                    266:      $                   Z( LDZ, * )
                    267: *     ..
                    268: *
                    269: *  =====================================================================
                    270: *
                    271: *     .. Parameters ..
                    272:       COMPLEX*16         CONE, CZERO
                    273:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ),
                    274:      $                   CZERO = ( 0.0D+0, 0.0D+0 ) )
                    275: *     ..
                    276: *     .. Local Scalars ..
                    277:       LOGICAL            LQUERY, UPPER, WANTZ
                    278:       CHARACTER          VECT
                    279:       INTEGER            IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLRWK,
                    280:      $                   LLWK2, LRWMIN, LWMIN
                    281: *     ..
                    282: *     .. External Functions ..
                    283:       LOGICAL            LSAME
                    284:       EXTERNAL           LSAME
                    285: *     ..
                    286: *     .. External Subroutines ..
                    287:       EXTERNAL           DSTERF, XERBLA, ZGEMM, ZHBGST, ZHBTRD, ZLACPY,
                    288:      $                   ZPBSTF, ZSTEDC
                    289: *     ..
                    290: *     .. Executable Statements ..
                    291: *
                    292: *     Test the input parameters.
                    293: *
                    294:       WANTZ = LSAME( JOBZ, 'V' )
                    295:       UPPER = LSAME( UPLO, 'U' )
                    296:       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
                    297: *
                    298:       INFO = 0
                    299:       IF( N.LE.1 ) THEN
1.8       bertrand  300:          LWMIN = 1+N
                    301:          LRWMIN = 1+N
1.1       bertrand  302:          LIWMIN = 1
                    303:       ELSE IF( WANTZ ) THEN
                    304:          LWMIN = 2*N**2
                    305:          LRWMIN = 1 + 5*N + 2*N**2
                    306:          LIWMIN = 3 + 5*N
                    307:       ELSE
                    308:          LWMIN = N
                    309:          LRWMIN = N
                    310:          LIWMIN = 1
                    311:       END IF
                    312:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
                    313:          INFO = -1
                    314:       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
                    315:          INFO = -2
                    316:       ELSE IF( N.LT.0 ) THEN
                    317:          INFO = -3
                    318:       ELSE IF( KA.LT.0 ) THEN
                    319:          INFO = -4
                    320:       ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
                    321:          INFO = -5
                    322:       ELSE IF( LDAB.LT.KA+1 ) THEN
                    323:          INFO = -7
                    324:       ELSE IF( LDBB.LT.KB+1 ) THEN
                    325:          INFO = -9
                    326:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
                    327:          INFO = -12
                    328:       END IF
                    329: *
                    330:       IF( INFO.EQ.0 ) THEN
                    331:          WORK( 1 ) = LWMIN
                    332:          RWORK( 1 ) = LRWMIN
                    333:          IWORK( 1 ) = LIWMIN
                    334: *
                    335:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
                    336:             INFO = -14
                    337:          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
                    338:             INFO = -16
                    339:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
                    340:             INFO = -18
                    341:          END IF
                    342:       END IF
                    343: *
                    344:       IF( INFO.NE.0 ) THEN
                    345:          CALL XERBLA( 'ZHBGVD', -INFO )
                    346:          RETURN
                    347:       ELSE IF( LQUERY ) THEN
                    348:          RETURN
                    349:       END IF
                    350: *
                    351: *     Quick return if possible
                    352: *
                    353:       IF( N.EQ.0 )
                    354:      $   RETURN
                    355: *
                    356: *     Form a split Cholesky factorization of B.
                    357: *
                    358:       CALL ZPBSTF( UPLO, N, KB, BB, LDBB, INFO )
                    359:       IF( INFO.NE.0 ) THEN
                    360:          INFO = N + INFO
                    361:          RETURN
                    362:       END IF
                    363: *
                    364: *     Transform problem to standard eigenvalue problem.
                    365: *
                    366:       INDE = 1
                    367:       INDWRK = INDE + N
                    368:       INDWK2 = 1 + N*N
                    369:       LLWK2 = LWORK - INDWK2 + 2
                    370:       LLRWK = LRWORK - INDWRK + 2
                    371:       CALL ZHBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ,
1.15      bertrand  372:      $             WORK, RWORK, IINFO )
1.1       bertrand  373: *
                    374: *     Reduce Hermitian band matrix to tridiagonal form.
                    375: *
                    376:       IF( WANTZ ) THEN
                    377:          VECT = 'U'
                    378:       ELSE
                    379:          VECT = 'N'
                    380:       END IF
                    381:       CALL ZHBTRD( VECT, UPLO, N, KA, AB, LDAB, W, RWORK( INDE ), Z,
                    382:      $             LDZ, WORK, IINFO )
                    383: *
                    384: *     For eigenvalues only, call DSTERF.  For eigenvectors, call ZSTEDC.
                    385: *
                    386:       IF( .NOT.WANTZ ) THEN
                    387:          CALL DSTERF( N, W, RWORK( INDE ), INFO )
                    388:       ELSE
                    389:          CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
                    390:      $                LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
                    391:      $                INFO )
                    392:          CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
                    393:      $               WORK( INDWK2 ), N )
                    394:          CALL ZLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
                    395:       END IF
                    396: *
                    397:       WORK( 1 ) = LWMIN
                    398:       RWORK( 1 ) = LRWMIN
                    399:       IWORK( 1 ) = LIWMIN
                    400:       RETURN
                    401: *
                    402: *     End of ZHBGVD
                    403: *
                    404:       END

CVSweb interface <joel.bertrand@systella.fr>