Return to dspgvd.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence Lapack.

    1: *> \brief \b DSPGVD
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DSPGVD + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dspgvd.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dspgvd.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dspgvd.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DSPGVD( ITYPE, JOBZ, UPLO, N, AP, BP, W, Z, LDZ, WORK,
   22: *                          LWORK, IWORK, LIWORK, INFO )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          JOBZ, UPLO
   26: *       INTEGER            INFO, ITYPE, LDZ, LIWORK, LWORK, N
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       INTEGER            IWORK( * )
   30: *       DOUBLE PRECISION   AP( * ), BP( * ), W( * ), WORK( * ),
   31: *      $                   Z( LDZ, * )
   32: *       ..
   33: *  
   34: *
   35: *> \par Purpose:
   36: *  =============
   37: *>
   38: *> \verbatim
   39: *>
   40: *> DSPGVD computes all the eigenvalues, and optionally, the eigenvectors
   41: *> of a real generalized symmetric-definite eigenproblem, of the form
   42: *> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
   43: *> B are assumed to be symmetric, stored in packed format, and B is also
   44: *> positive definite.
   45: *> If eigenvectors are desired, it uses a divide and conquer algorithm.
   46: *>
   47: *> The divide and conquer algorithm makes very mild assumptions about
   48: *> floating point arithmetic. It will work on machines with a guard
   49: *> digit in add/subtract, or on those binary machines without guard
   50: *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
   51: *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
   52: *> without guard digits, but we know of none.
   53: *> \endverbatim
   54: *
   55: *  Arguments:
   56: *  ==========
   57: *
   58: *> \param[in] ITYPE
   59: *> \verbatim
   60: *>          ITYPE is INTEGER
   61: *>          Specifies the problem type to be solved:
   62: *>          = 1:  A*x = (lambda)*B*x
   63: *>          = 2:  A*B*x = (lambda)*x
   64: *>          = 3:  B*A*x = (lambda)*x
   65: *> \endverbatim
   66: *>
   67: *> \param[in] JOBZ
   68: *> \verbatim
   69: *>          JOBZ is CHARACTER*1
   70: *>          = 'N':  Compute eigenvalues only;
   71: *>          = 'V':  Compute eigenvalues and eigenvectors.
   72: *> \endverbatim
   73: *>
   74: *> \param[in] UPLO
   75: *> \verbatim
   76: *>          UPLO is CHARACTER*1
   77: *>          = 'U':  Upper triangles of A and B are stored;
   78: *>          = 'L':  Lower triangles of A and B are stored.
   79: *> \endverbatim
   80: *>
   81: *> \param[in] N
   82: *> \verbatim
   83: *>          N is INTEGER
   84: *>          The order of the matrices A and B.  N >= 0.
   85: *> \endverbatim
   86: *>
   87: *> \param[in,out] AP
   88: *> \verbatim
   89: *>          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
   90: *>          On entry, the upper or lower triangle of the symmetric matrix
   91: *>          A, packed columnwise in a linear array.  The j-th column of A
   92: *>          is stored in the array AP as follows:
   93: *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
   94: *>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
   95: *>
   96: *>          On exit, the contents of AP are destroyed.
   97: *> \endverbatim
   98: *>
   99: *> \param[in,out] BP
  100: *> \verbatim
  101: *>          BP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
  102: *>          On entry, the upper or lower triangle of the symmetric matrix
  103: *>          B, packed columnwise in a linear array.  The j-th column of B
  104: *>          is stored in the array BP as follows:
  105: *>          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
  106: *>          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
  107: *>
  108: *>          On exit, the triangular factor U or L from the Cholesky
  109: *>          factorization B = U**T*U or B = L*L**T, in the same storage
  110: *>          format as B.
  111: *> \endverbatim
  112: *>
  113: *> \param[out] W
  114: *> \verbatim
  115: *>          W is DOUBLE PRECISION array, dimension (N)
  116: *>          If INFO = 0, the eigenvalues in ascending order.
  117: *> \endverbatim
  118: *>
  119: *> \param[out] Z
  120: *> \verbatim
  121: *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
  122: *>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
  123: *>          eigenvectors.  The eigenvectors are normalized as follows:
  124: *>          if ITYPE = 1 or 2, Z**T*B*Z = I;
  125: *>          if ITYPE = 3, Z**T*inv(B)*Z = I.
  126: *>          If JOBZ = 'N', then Z is not referenced.
  127: *> \endverbatim
  128: *>
  129: *> \param[in] LDZ
  130: *> \verbatim
  131: *>          LDZ is INTEGER
  132: *>          The leading dimension of the array Z.  LDZ >= 1, and if
  133: *>          JOBZ = 'V', LDZ >= max(1,N).
  134: *> \endverbatim
  135: *>
  136: *> \param[out] WORK
  137: *> \verbatim
  138: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  139: *>          On exit, if INFO = 0, WORK(1) returns the required LWORK.
  140: *> \endverbatim
  141: *>
  142: *> \param[in] LWORK
  143: *> \verbatim
  144: *>          LWORK is INTEGER
  145: *>          The dimension of the array WORK.
  146: *>          If N <= 1,               LWORK >= 1.
  147: *>          If JOBZ = 'N' and N > 1, LWORK >= 2*N.
  148: *>          If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
  149: *>
  150: *>          If LWORK = -1, then a workspace query is assumed; the routine
  151: *>          only calculates the required sizes of the WORK and IWORK
  152: *>          arrays, returns these values as the first entries of the WORK
  153: *>          and IWORK arrays, and no error message related to LWORK or
  154: *>          LIWORK is issued by XERBLA.
  155: *> \endverbatim
  156: *>
  157: *> \param[out] IWORK
  158: *> \verbatim
  159: *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
  160: *>          On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
  161: *> \endverbatim
  162: *>
  163: *> \param[in] LIWORK
  164: *> \verbatim
  165: *>          LIWORK is INTEGER
  166: *>          The dimension of the array IWORK.
  167: *>          If JOBZ  = 'N' or N <= 1, LIWORK >= 1.
  168: *>          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
  169: *>
  170: *>          If LIWORK = -1, then a workspace query is assumed; the
  171: *>          routine only calculates the required sizes of the WORK and
  172: *>          IWORK arrays, returns these values as the first entries of
  173: *>          the WORK and IWORK arrays, and no error message related to
  174: *>          LWORK or LIWORK is issued by XERBLA.
  175: *> \endverbatim
  176: *>
  177: *> \param[out] INFO
  178: *> \verbatim
  179: *>          INFO is INTEGER
  180: *>          = 0:  successful exit
  181: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  182: *>          > 0:  DPPTRF or DSPEVD returned an error code:
  183: *>             <= N:  if INFO = i, DSPEVD failed to converge;
  184: *>                    i off-diagonal elements of an intermediate
  185: *>                    tridiagonal form did not converge to zero;
  186: *>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
  187: *>                    minor of order i of B is not positive definite.
  188: *>                    The factorization of B could not be completed and
  189: *>                    no eigenvalues or eigenvectors were computed.
  190: *> \endverbatim
  191: *
  192: *  Authors:
  193: *  ========
  194: *
  195: *> \author Univ. of Tennessee 
  196: *> \author Univ. of California Berkeley 
  197: *> \author Univ. of Colorado Denver 
  198: *> \author NAG Ltd. 
  199: *
  200: *> \date November 2015
  201: *
  202: *> \ingroup doubleOTHEReigen
  203: *
  204: *> \par Contributors:
  205: *  ==================
  206: *>
  207: *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
  208: *
  209: *  =====================================================================
  210:       SUBROUTINE DSPGVD( ITYPE, JOBZ, UPLO, N, AP, BP, W, Z, LDZ, WORK,
  211:      $                   LWORK, IWORK, LIWORK, INFO )
  212: *
  213: *  -- LAPACK driver routine (version 3.6.0) --
  214: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  215: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  216: *     November 2015
  217: *
  218: *     .. Scalar Arguments ..
  219:       CHARACTER          JOBZ, UPLO
  220:       INTEGER            INFO, ITYPE, LDZ, LIWORK, LWORK, N
  221: *     ..
  222: *     .. Array Arguments ..
  223:       INTEGER            IWORK( * )
  224:       DOUBLE PRECISION   AP( * ), BP( * ), W( * ), WORK( * ),
  225:      $                   Z( LDZ, * )
  226: *     ..
  227: *
  228: *  =====================================================================
  229: *
  230: *     .. Local Scalars ..
  231:       LOGICAL            LQUERY, UPPER, WANTZ
  232:       CHARACTER          TRANS
  233:       INTEGER            J, LIWMIN, LWMIN, NEIG
  234: *     ..
  235: *     .. External Functions ..
  236:       LOGICAL            LSAME
  237:       EXTERNAL           LSAME
  238: *     ..
  239: *     .. External Subroutines ..
  240:       EXTERNAL           DPPTRF, DSPEVD, DSPGST, DTPMV, DTPSV, XERBLA
  241: *     ..
  242: *     .. Intrinsic Functions ..
  243:       INTRINSIC          DBLE, MAX
  244: *     ..
  245: *     .. Executable Statements ..
  246: *
  247: *     Test the input parameters.
  248: *
  249:       WANTZ = LSAME( JOBZ, 'V' )
  250:       UPPER = LSAME( UPLO, 'U' )
  251:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
  252: *
  253:       INFO = 0
  254:       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
  255:          INFO = -1
  256:       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
  257:          INFO = -2
  258:       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
  259:          INFO = -3
  260:       ELSE IF( N.LT.0 ) THEN
  261:          INFO = -4
  262:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
  263:          INFO = -9
  264:       END IF
  265: *
  266:       IF( INFO.EQ.0 ) THEN
  267:          IF( N.LE.1 ) THEN
  268:             LIWMIN = 1
  269:             LWMIN = 1
  270:          ELSE
  271:             IF( WANTZ ) THEN
  272:                LIWMIN = 3 + 5*N
  273:                LWMIN = 1 + 6*N + 2*N**2
  274:             ELSE
  275:                LIWMIN = 1
  276:                LWMIN = 2*N
  277:             END IF
  278:          END IF
  279:          WORK( 1 ) = LWMIN
  280:          IWORK( 1 ) = LIWMIN
  281:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  282:             INFO = -11
  283:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
  284:             INFO = -13
  285:          END IF
  286:       END IF
  287: *
  288:       IF( INFO.NE.0 ) THEN
  289:          CALL XERBLA( 'DSPGVD', -INFO )
  290:          RETURN
  291:       ELSE IF( LQUERY ) THEN
  292:          RETURN
  293:       END IF
  294: *
  295: *     Quick return if possible
  296: *
  297:       IF( N.EQ.0 )
  298:      $   RETURN
  299: *
  300: *     Form a Cholesky factorization of BP.
  301: *
  302:       CALL DPPTRF( UPLO, N, BP, INFO )
  303:       IF( INFO.NE.0 ) THEN
  304:          INFO = N + INFO
  305:          RETURN
  306:       END IF
  307: *
  308: *     Transform problem to standard eigenvalue problem and solve.
  309: *
  310:       CALL DSPGST( ITYPE, UPLO, N, AP, BP, INFO )
  311:       CALL DSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, IWORK,
  312:      $             LIWORK, INFO )
  313:       LWMIN = MAX( DBLE( LWMIN ), DBLE( WORK( 1 ) ) )
  314:       LIWMIN = MAX( DBLE( LIWMIN ), DBLE( IWORK( 1 ) ) )
  315: *
  316:       IF( WANTZ ) THEN
  317: *
  318: *        Backtransform eigenvectors to the original problem.
  319: *
  320:          NEIG = N
  321:          IF( INFO.GT.0 )
  322:      $      NEIG = INFO - 1
  323:          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
  324: *
  325: *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
  326: *           backtransform eigenvectors: x = inv(L)**T *y or inv(U)*y
  327: *
  328:             IF( UPPER ) THEN
  329:                TRANS = 'N'
  330:             ELSE
  331:                TRANS = 'T'
  332:             END IF
  333: *
  334:             DO 10 J = 1, NEIG
  335:                CALL DTPSV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ),
  336:      $                     1 )
  337:    10       CONTINUE
  338: *
  339:          ELSE IF( ITYPE.EQ.3 ) THEN
  340: *
  341: *           For B*A*x=(lambda)*x;
  342: *           backtransform eigenvectors: x = L*y or U**T *y
  343: *
  344:             IF( UPPER ) THEN
  345:                TRANS = 'T'
  346:             ELSE
  347:                TRANS = 'N'
  348:             END IF
  349: *
  350:             DO 20 J = 1, NEIG
  351:                CALL DTPMV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ),
  352:      $                     1 )
  353:    20       CONTINUE
  354:          END IF
  355:       END IF
  356: *
  357:       WORK( 1 ) = LWMIN
  358:       IWORK( 1 ) = LIWMIN
  359: *
  360:       RETURN
  361: *
  362: *     End of DSPGVD
  363: *
  364:       END