Return to dsbgvd.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

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