Return to zhbgvd.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZHBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
    2:      $                   Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK,
    3:      $                   LIWORK, INFO )
    4: *
    5: *  -- LAPACK driver routine (version 3.2) --
    6: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    7: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    8: *     November 2006
    9: *
   10: *     .. Scalar Arguments ..
   11:       CHARACTER          JOBZ, UPLO
   12:       INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LRWORK,
   13:      $                   LWORK, N
   14: *     ..
   15: *     .. Array Arguments ..
   16:       INTEGER            IWORK( * )
   17:       DOUBLE PRECISION   RWORK( * ), W( * )
   18:       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
   19:      $                   Z( LDZ, * )
   20: *     ..
   21: *
   22: *  Purpose
   23: *  =======
   24: *
   25: *  ZHBGVD computes all the eigenvalues, and optionally, the eigenvectors
   26: *  of a complex generalized Hermitian-definite banded eigenproblem, of
   27: *  the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian
   28: *  and banded, and B is also positive definite.  If eigenvectors are
   29: *  desired, it uses a divide and conquer algorithm.
   30: *
   31: *  The divide and conquer algorithm makes very mild assumptions about
   32: *  floating point arithmetic. It will work on machines with a guard
   33: *  digit in add/subtract, or on those binary machines without guard
   34: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
   35: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
   36: *  without guard digits, but we know of none.
   37: *
   38: *  Arguments
   39: *  =========
   40: *
   41: *  JOBZ    (input) CHARACTER*1
   42: *          = 'N':  Compute eigenvalues only;
   43: *          = 'V':  Compute eigenvalues and eigenvectors.
   44: *
   45: *  UPLO    (input) CHARACTER*1
   46: *          = 'U':  Upper triangles of A and B are stored;
   47: *          = 'L':  Lower triangles of A and B are stored.
   48: *
   49: *  N       (input) INTEGER
   50: *          The order of the matrices A and B.  N >= 0.
   51: *
   52: *  KA      (input) INTEGER
   53: *          The number of superdiagonals of the matrix A if UPLO = 'U',
   54: *          or the number of subdiagonals if UPLO = 'L'. KA >= 0.
   55: *
   56: *  KB      (input) INTEGER
   57: *          The number of superdiagonals of the matrix B if UPLO = 'U',
   58: *          or the number of subdiagonals if UPLO = 'L'. KB >= 0.
   59: *
   60: *  AB      (input/output) COMPLEX*16 array, dimension (LDAB, N)
   61: *          On entry, the upper or lower triangle of the Hermitian band
   62: *          matrix A, stored in the first ka+1 rows of the array.  The
   63: *          j-th column of A is stored in the j-th column of the array AB
   64: *          as follows:
   65: *          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
   66: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
   67: *
   68: *          On exit, the contents of AB are destroyed.
   69: *
   70: *  LDAB    (input) INTEGER
   71: *          The leading dimension of the array AB.  LDAB >= KA+1.
   72: *
   73: *  BB      (input/output) COMPLEX*16 array, dimension (LDBB, N)
   74: *          On entry, the upper or lower triangle of the Hermitian band
   75: *          matrix B, stored in the first kb+1 rows of the array.  The
   76: *          j-th column of B is stored in the j-th column of the array BB
   77: *          as follows:
   78: *          if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
   79: *          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
   80: *
   81: *          On exit, the factor S from the split Cholesky factorization
   82: *          B = S**H*S, as returned by ZPBSTF.
   83: *
   84: *  LDBB    (input) INTEGER
   85: *          The leading dimension of the array BB.  LDBB >= KB+1.
   86: *
   87: *  W       (output) DOUBLE PRECISION array, dimension (N)
   88: *          If INFO = 0, the eigenvalues in ascending order.
   89: *
   90: *  Z       (output) COMPLEX*16 array, dimension (LDZ, N)
   91: *          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
   92: *          eigenvectors, with the i-th column of Z holding the
   93: *          eigenvector associated with W(i). The eigenvectors are
   94: *          normalized so that Z**H*B*Z = I.
   95: *          If JOBZ = 'N', then Z is not referenced.
   96: *
   97: *  LDZ     (input) INTEGER
   98: *          The leading dimension of the array Z.  LDZ >= 1, and if
   99: *          JOBZ = 'V', LDZ >= N.
  100: *
  101: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
  102: *          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
  103: *
  104: *  LWORK   (input) INTEGER
  105: *          The dimension of the array WORK.
  106: *          If N <= 1,               LWORK >= 1.
  107: *          If JOBZ = 'N' and N > 1, LWORK >= N.
  108: *          If JOBZ = 'V' and N > 1, LWORK >= 2*N**2.
  109: *
  110: *          If LWORK = -1, then a workspace query is assumed; the routine
  111: *          only calculates the optimal sizes of the WORK, RWORK and
  112: *          IWORK arrays, returns these values as the first entries of
  113: *          the WORK, RWORK and IWORK arrays, and no error message
  114: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
  115: *
  116: *  RWORK   (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
  117: *          On exit, if INFO=0, RWORK(1) returns the optimal LRWORK.
  118: *
  119: *  LRWORK  (input) INTEGER
  120: *          The dimension of array RWORK.
  121: *          If N <= 1,               LRWORK >= 1.
  122: *          If JOBZ = 'N' and N > 1, LRWORK >= N.
  123: *          If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
  124: *
  125: *          If LRWORK = -1, then a workspace query is assumed; the
  126: *          routine only calculates the optimal sizes of the WORK, RWORK
  127: *          and IWORK arrays, returns these values as the first entries
  128: *          of the WORK, RWORK and IWORK arrays, and no error message
  129: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
  130: *
  131: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
  132: *          On exit, if INFO=0, IWORK(1) returns the optimal LIWORK.
  133: *
  134: *  LIWORK  (input) INTEGER
  135: *          The dimension of array IWORK.
  136: *          If JOBZ = 'N' or N <= 1, LIWORK >= 1.
  137: *          If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
  138: *
  139: *          If LIWORK = -1, then a workspace query is assumed; the
  140: *          routine only calculates the optimal sizes of the WORK, RWORK
  141: *          and IWORK arrays, returns these values as the first entries
  142: *          of the WORK, RWORK and IWORK arrays, and no error message
  143: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
  144: *
  145: *  INFO    (output) INTEGER
  146: *          = 0:  successful exit
  147: *          < 0:  if INFO = -i, the i-th argument had an illegal value
  148: *          > 0:  if INFO = i, and i is:
  149: *             <= N:  the algorithm failed to converge:
  150: *                    i off-diagonal elements of an intermediate
  151: *                    tridiagonal form did not converge to zero;
  152: *             > N:   if INFO = N + i, for 1 <= i <= N, then ZPBSTF
  153: *                    returned INFO = i: B is not positive definite.
  154: *                    The factorization of B could not be completed and
  155: *                    no eigenvalues or eigenvectors were computed.
  156: *
  157: *  Further Details
  158: *  ===============
  159: *
  160: *  Based on contributions by
  161: *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
  162: *
  163: *  =====================================================================
  164: *
  165: *     .. Parameters ..
  166:       COMPLEX*16         CONE, CZERO
  167:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ),
  168:      $                   CZERO = ( 0.0D+0, 0.0D+0 ) )
  169: *     ..
  170: *     .. Local Scalars ..
  171:       LOGICAL            LQUERY, UPPER, WANTZ
  172:       CHARACTER          VECT
  173:       INTEGER            IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLRWK,
  174:      $                   LLWK2, LRWMIN, LWMIN
  175: *     ..
  176: *     .. External Functions ..
  177:       LOGICAL            LSAME
  178:       EXTERNAL           LSAME
  179: *     ..
  180: *     .. External Subroutines ..
  181:       EXTERNAL           DSTERF, XERBLA, ZGEMM, ZHBGST, ZHBTRD, ZLACPY,
  182:      $                   ZPBSTF, ZSTEDC
  183: *     ..
  184: *     .. Executable Statements ..
  185: *
  186: *     Test the input parameters.
  187: *
  188:       WANTZ = LSAME( JOBZ, 'V' )
  189:       UPPER = LSAME( UPLO, 'U' )
  190:       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
  191: *
  192:       INFO = 0
  193:       IF( N.LE.1 ) THEN
  194:          LWMIN = 1
  195:          LRWMIN = 1
  196:          LIWMIN = 1
  197:       ELSE IF( WANTZ ) THEN
  198:          LWMIN = 2*N**2
  199:          LRWMIN = 1 + 5*N + 2*N**2
  200:          LIWMIN = 3 + 5*N
  201:       ELSE
  202:          LWMIN = N
  203:          LRWMIN = N
  204:          LIWMIN = 1
  205:       END IF
  206:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
  207:          INFO = -1
  208:       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
  209:          INFO = -2
  210:       ELSE IF( N.LT.0 ) THEN
  211:          INFO = -3
  212:       ELSE IF( KA.LT.0 ) THEN
  213:          INFO = -4
  214:       ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
  215:          INFO = -5
  216:       ELSE IF( LDAB.LT.KA+1 ) THEN
  217:          INFO = -7
  218:       ELSE IF( LDBB.LT.KB+1 ) THEN
  219:          INFO = -9
  220:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
  221:          INFO = -12
  222:       END IF
  223: *
  224:       IF( INFO.EQ.0 ) THEN
  225:          WORK( 1 ) = LWMIN
  226:          RWORK( 1 ) = LRWMIN
  227:          IWORK( 1 ) = LIWMIN
  228: *
  229:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  230:             INFO = -14
  231:          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
  232:             INFO = -16
  233:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
  234:             INFO = -18
  235:          END IF
  236:       END IF
  237: *
  238:       IF( INFO.NE.0 ) THEN
  239:          CALL XERBLA( 'ZHBGVD', -INFO )
  240:          RETURN
  241:       ELSE IF( LQUERY ) THEN
  242:          RETURN
  243:       END IF
  244: *
  245: *     Quick return if possible
  246: *
  247:       IF( N.EQ.0 )
  248:      $   RETURN
  249: *
  250: *     Form a split Cholesky factorization of B.
  251: *
  252:       CALL ZPBSTF( UPLO, N, KB, BB, LDBB, INFO )
  253:       IF( INFO.NE.0 ) THEN
  254:          INFO = N + INFO
  255:          RETURN
  256:       END IF
  257: *
  258: *     Transform problem to standard eigenvalue problem.
  259: *
  260:       INDE = 1
  261:       INDWRK = INDE + N
  262:       INDWK2 = 1 + N*N
  263:       LLWK2 = LWORK - INDWK2 + 2
  264:       LLRWK = LRWORK - INDWRK + 2
  265:       CALL ZHBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ,
  266:      $             WORK, RWORK( INDWRK ), IINFO )
  267: *
  268: *     Reduce Hermitian band matrix to tridiagonal form.
  269: *
  270:       IF( WANTZ ) THEN
  271:          VECT = 'U'
  272:       ELSE
  273:          VECT = 'N'
  274:       END IF
  275:       CALL ZHBTRD( VECT, UPLO, N, KA, AB, LDAB, W, RWORK( INDE ), Z,
  276:      $             LDZ, WORK, IINFO )
  277: *
  278: *     For eigenvalues only, call DSTERF.  For eigenvectors, call ZSTEDC.
  279: *
  280:       IF( .NOT.WANTZ ) THEN
  281:          CALL DSTERF( N, W, RWORK( INDE ), INFO )
  282:       ELSE
  283:          CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
  284:      $                LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
  285:      $                INFO )
  286:          CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
  287:      $               WORK( INDWK2 ), N )
  288:          CALL ZLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
  289:       END IF
  290: *
  291:       WORK( 1 ) = LWMIN
  292:       RWORK( 1 ) = LRWMIN
  293:       IWORK( 1 ) = LIWMIN
  294:       RETURN
  295: *
  296: *     End of ZHBGVD
  297: *
  298:       END