Return to dspgvd.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

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