File:  [local] / rpl / lapack / lapack / zhbgvd.f
Revision 1.16: download - view: text, annotated - select for diffs - revision graph
Sat Aug 27 15:34:49 2016 UTC (7 years, 8 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_25, HEAD
Cohérence Lapack.

    1: *> \brief \b ZHBGVD
    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 June 2016
  242: *
  243: *> \ingroup complex16OTHEReigen
  244: *
  245: *> \par Contributors:
  246: *  ==================
  247: *>
  248: *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
  249: *
  250: *  =====================================================================
  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: *
  255: *  -- LAPACK driver routine (version 3.6.1) --
  256: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  257: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  258: *     June 2016
  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
  303:          LWMIN = 1+N
  304:          LRWMIN = 1+N
  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, 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>