Return to dsygvd.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
    2:      $                   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, ITYPE, LDA, LDB, LIWORK, LWORK, N
   12: *     ..
   13: *     .. Array Arguments ..
   14:       INTEGER            IWORK( * )
   15:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  DSYGVD computes all the eigenvalues, and optionally, the eigenvectors
   22: *  of a real generalized symmetric-definite eigenproblem, of the form
   23: *  A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
   24: *  B are assumed to be symmetric and B is also positive definite.
   25: *  If eigenvectors are desired, it uses a divide and conquer algorithm.
   26: *
   27: *  The divide and conquer algorithm makes very mild assumptions about
   28: *  floating point arithmetic. It will work on machines with a guard
   29: *  digit in add/subtract, or on those binary machines without guard
   30: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
   31: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
   32: *  without guard digits, but we know of none.
   33: *
   34: *  Arguments
   35: *  =========
   36: *
   37: *  ITYPE   (input) INTEGER
   38: *          Specifies the problem type to be solved:
   39: *          = 1:  A*x = (lambda)*B*x
   40: *          = 2:  A*B*x = (lambda)*x
   41: *          = 3:  B*A*x = (lambda)*x
   42: *
   43: *  JOBZ    (input) CHARACTER*1
   44: *          = 'N':  Compute eigenvalues only;
   45: *          = 'V':  Compute eigenvalues and eigenvectors.
   46: *
   47: *  UPLO    (input) CHARACTER*1
   48: *          = 'U':  Upper triangles of A and B are stored;
   49: *          = 'L':  Lower triangles of A and B are stored.
   50: *
   51: *  N       (input) INTEGER
   52: *          The order of the matrices A and B.  N >= 0.
   53: *
   54: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
   55: *          On entry, the symmetric matrix A.  If UPLO = 'U', the
   56: *          leading N-by-N upper triangular part of A contains the
   57: *          upper triangular part of the matrix A.  If UPLO = 'L',
   58: *          the leading N-by-N lower triangular part of A contains
   59: *          the lower triangular part of the matrix A.
   60: *
   61: *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
   62: *          matrix Z of eigenvectors.  The eigenvectors are normalized
   63: *          as follows:
   64: *          if ITYPE = 1 or 2, Z**T*B*Z = I;
   65: *          if ITYPE = 3, Z**T*inv(B)*Z = I.
   66: *          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
   67: *          or the lower triangle (if UPLO='L') of A, including the
   68: *          diagonal, is destroyed.
   69: *
   70: *  LDA     (input) INTEGER
   71: *          The leading dimension of the array A.  LDA >= max(1,N).
   72: *
   73: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
   74: *          On entry, the symmetric matrix B.  If UPLO = 'U', the
   75: *          leading N-by-N upper triangular part of B contains the
   76: *          upper triangular part of the matrix B.  If UPLO = 'L',
   77: *          the leading N-by-N lower triangular part of B contains
   78: *          the lower triangular part of the matrix B.
   79: *
   80: *          On exit, if INFO <= N, the part of B containing the matrix is
   81: *          overwritten by the triangular factor U or L from the Cholesky
   82: *          factorization B = U**T*U or B = L*L**T.
   83: *
   84: *  LDB     (input) INTEGER
   85: *          The leading dimension of the array B.  LDB >= max(1,N).
   86: *
   87: *  W       (output) DOUBLE PRECISION array, dimension (N)
   88: *          If INFO = 0, the eigenvalues in ascending order.
   89: *
   90: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
   91: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   92: *
   93: *  LWORK   (input) INTEGER
   94: *          The dimension of the array WORK.
   95: *          If N <= 1,               LWORK >= 1.
   96: *          If JOBZ = 'N' and N > 1, LWORK >= 2*N+1.
   97: *          If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
   98: *
   99: *          If LWORK = -1, then a workspace query is assumed; the routine
  100: *          only calculates the optimal sizes of the WORK and IWORK
  101: *          arrays, returns these values as the first entries of the WORK
  102: *          and IWORK arrays, and no error message related to LWORK or
  103: *          LIWORK is issued by XERBLA.
  104: *
  105: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
  106: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
  107: *
  108: *  LIWORK  (input) INTEGER
  109: *          The dimension of the array IWORK.
  110: *          If N <= 1,                LIWORK >= 1.
  111: *          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
  112: *          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
  113: *
  114: *          If LIWORK = -1, then a workspace query is assumed; the
  115: *          routine only calculates the optimal sizes of the WORK and
  116: *          IWORK arrays, returns these values as the first entries of
  117: *          the WORK and IWORK arrays, and no error message related to
  118: *          LWORK or LIWORK is issued by XERBLA.
  119: *
  120: *  INFO    (output) INTEGER
  121: *          = 0:  successful exit
  122: *          < 0:  if INFO = -i, the i-th argument had an illegal value
  123: *          > 0:  DPOTRF or DSYEVD returned an error code:
  124: *             <= N:  if INFO = i and JOBZ = 'N', then the algorithm
  125: *                    failed to converge; i off-diagonal elements of an
  126: *                    intermediate tridiagonal form did not converge to
  127: *                    zero;
  128: *                    if INFO = i and JOBZ = 'V', then the algorithm
  129: *                    failed to compute an eigenvalue while working on
  130: *                    the submatrix lying in rows and columns INFO/(N+1)
  131: *                    through mod(INFO,N+1);
  132: *             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
  133: *                    minor of order i of B is not positive definite.
  134: *                    The factorization of B could not be completed and
  135: *                    no eigenvalues or eigenvectors were computed.
  136: *
  137: *  Further Details
  138: *  ===============
  139: *
  140: *  Based on contributions by
  141: *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
  142: *
  143: *  Modified so that no backsubstitution is performed if DSYEVD fails to
  144: *  converge (NEIG in old code could be greater than N causing out of
  145: *  bounds reference to A - reported by Ralf Meyer).  Also corrected the
  146: *  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
  147: *  =====================================================================
  148: *
  149: *     .. Parameters ..
  150:       DOUBLE PRECISION   ONE
  151:       PARAMETER          ( ONE = 1.0D+0 )
  152: *     ..
  153: *     .. Local Scalars ..
  154:       LOGICAL            LQUERY, UPPER, WANTZ
  155:       CHARACTER          TRANS
  156:       INTEGER            LIOPT, LIWMIN, LOPT, LWMIN
  157: *     ..
  158: *     .. External Functions ..
  159:       LOGICAL            LSAME
  160:       EXTERNAL           LSAME
  161: *     ..
  162: *     .. External Subroutines ..
  163:       EXTERNAL           DPOTRF, DSYEVD, DSYGST, DTRMM, DTRSM, XERBLA
  164: *     ..
  165: *     .. Intrinsic Functions ..
  166:       INTRINSIC          DBLE, MAX
  167: *     ..
  168: *     .. Executable Statements ..
  169: *
  170: *     Test the input parameters.
  171: *
  172:       WANTZ = LSAME( JOBZ, 'V' )
  173:       UPPER = LSAME( UPLO, 'U' )
  174:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
  175: *
  176:       INFO = 0
  177:       IF( N.LE.1 ) THEN
  178:          LIWMIN = 1
  179:          LWMIN = 1
  180:       ELSE IF( WANTZ ) THEN
  181:          LIWMIN = 3 + 5*N
  182:          LWMIN = 1 + 6*N + 2*N**2
  183:       ELSE
  184:          LIWMIN = 1
  185:          LWMIN = 2*N + 1
  186:       END IF
  187:       LOPT = LWMIN
  188:       LIOPT = LIWMIN
  189:       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
  190:          INFO = -1
  191:       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
  192:          INFO = -2
  193:       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
  194:          INFO = -3
  195:       ELSE IF( N.LT.0 ) THEN
  196:          INFO = -4
  197:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  198:          INFO = -6
  199:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  200:          INFO = -8
  201:       END IF
  202: *
  203:       IF( INFO.EQ.0 ) THEN
  204:          WORK( 1 ) = LOPT
  205:          IWORK( 1 ) = LIOPT
  206: *
  207:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  208:             INFO = -11
  209:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
  210:             INFO = -13
  211:          END IF
  212:       END IF
  213: *
  214:       IF( INFO.NE.0 ) THEN
  215:          CALL XERBLA( 'DSYGVD', -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 Cholesky factorization of B.
  227: *
  228:       CALL DPOTRF( UPLO, N, B, LDB, INFO )
  229:       IF( INFO.NE.0 ) THEN
  230:          INFO = N + INFO
  231:          RETURN
  232:       END IF
  233: *
  234: *     Transform problem to standard eigenvalue problem and solve.
  235: *
  236:       CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
  237:       CALL DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK,
  238:      $             INFO )
  239:       LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) )
  240:       LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) )
  241: *
  242:       IF( WANTZ .AND. INFO.EQ.0 ) THEN
  243: *
  244: *        Backtransform eigenvectors to the original problem.
  245: *
  246:          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
  247: *
  248: *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
  249: *           backtransform eigenvectors: x = inv(L)'*y or inv(U)*y
  250: *
  251:             IF( UPPER ) THEN
  252:                TRANS = 'N'
  253:             ELSE
  254:                TRANS = 'T'
  255:             END IF
  256: *
  257:             CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
  258:      $                  B, LDB, A, LDA )
  259: *
  260:          ELSE IF( ITYPE.EQ.3 ) THEN
  261: *
  262: *           For B*A*x=(lambda)*x;
  263: *           backtransform eigenvectors: x = L*y or U'*y
  264: *
  265:             IF( UPPER ) THEN
  266:                TRANS = 'T'
  267:             ELSE
  268:                TRANS = 'N'
  269:             END IF
  270: *
  271:             CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
  272:      $                  B, LDB, A, LDA )
  273:          END IF
  274:       END IF
  275: *
  276:       WORK( 1 ) = LOPT
  277:       IWORK( 1 ) = LIOPT
  278: *
  279:       RETURN
  280: *
  281: *     End of DSYGVD
  282: *
  283:       END