Return to dpteqr.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

    1:       SUBROUTINE DPTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       CHARACTER          COMPZ
   10:       INTEGER            INFO, LDZ, N
   11: *     ..
   12: *     .. Array Arguments ..
   13:       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  DPTEQR computes all eigenvalues and, optionally, eigenvectors of a
   20: *  symmetric positive definite tridiagonal matrix by first factoring the
   21: *  matrix using DPTTRF, and then calling DBDSQR to compute the singular
   22: *  values of the bidiagonal factor.
   23: *
   24: *  This routine computes the eigenvalues of the positive definite
   25: *  tridiagonal matrix to high relative accuracy.  This means that if the
   26: *  eigenvalues range over many orders of magnitude in size, then the
   27: *  small eigenvalues and corresponding eigenvectors will be computed
   28: *  more accurately than, for example, with the standard QR method.
   29: *
   30: *  The eigenvectors of a full or band symmetric positive definite matrix
   31: *  can also be found if DSYTRD, DSPTRD, or DSBTRD has been used to
   32: *  reduce this matrix to tridiagonal form. (The reduction to tridiagonal
   33: *  form, however, may preclude the possibility of obtaining high
   34: *  relative accuracy in the small eigenvalues of the original matrix, if
   35: *  these eigenvalues range over many orders of magnitude.)
   36: *
   37: *  Arguments
   38: *  =========
   39: *
   40: *  COMPZ   (input) CHARACTER*1
   41: *          = 'N':  Compute eigenvalues only.
   42: *          = 'V':  Compute eigenvectors of original symmetric
   43: *                  matrix also.  Array Z contains the orthogonal
   44: *                  matrix used to reduce the original matrix to
   45: *                  tridiagonal form.
   46: *          = 'I':  Compute eigenvectors of tridiagonal matrix also.
   47: *
   48: *  N       (input) INTEGER
   49: *          The order of the matrix.  N >= 0.
   50: *
   51: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
   52: *          On entry, the n diagonal elements of the tridiagonal
   53: *          matrix.
   54: *          On normal exit, D contains the eigenvalues, in descending
   55: *          order.
   56: *
   57: *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
   58: *          On entry, the (n-1) subdiagonal elements of the tridiagonal
   59: *          matrix.
   60: *          On exit, E has been destroyed.
   61: *
   62: *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
   63: *          On entry, if COMPZ = 'V', the orthogonal matrix used in the
   64: *          reduction to tridiagonal form.
   65: *          On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
   66: *          original symmetric matrix;
   67: *          if COMPZ = 'I', the orthonormal eigenvectors of the
   68: *          tridiagonal matrix.
   69: *          If INFO > 0 on exit, Z contains the eigenvectors associated
   70: *          with only the stored eigenvalues.
   71: *          If  COMPZ = 'N', then Z is not referenced.
   72: *
   73: *  LDZ     (input) INTEGER
   74: *          The leading dimension of the array Z.  LDZ >= 1, and if
   75: *          COMPZ = 'V' or 'I', LDZ >= max(1,N).
   76: *
   77: *  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)
   78: *
   79: *  INFO    (output) INTEGER
   80: *          = 0:  successful exit.
   81: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
   82: *          > 0:  if INFO = i, and i is:
   83: *                <= N  the Cholesky factorization of the matrix could
   84: *                      not be performed because the i-th principal minor
   85: *                      was not positive definite.
   86: *                > N   the SVD algorithm failed to converge;
   87: *                      if INFO = N+i, i off-diagonal elements of the
   88: *                      bidiagonal factor did not converge to zero.
   89: *
   90: *  =====================================================================
   91: *
   92: *     .. Parameters ..
   93:       DOUBLE PRECISION   ZERO, ONE
   94:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
   95: *     ..
   96: *     .. External Functions ..
   97:       LOGICAL            LSAME
   98:       EXTERNAL           LSAME
   99: *     ..
  100: *     .. External Subroutines ..
  101:       EXTERNAL           DBDSQR, DLASET, DPTTRF, XERBLA
  102: *     ..
  103: *     .. Local Arrays ..
  104:       DOUBLE PRECISION   C( 1, 1 ), VT( 1, 1 )
  105: *     ..
  106: *     .. Local Scalars ..
  107:       INTEGER            I, ICOMPZ, NRU
  108: *     ..
  109: *     .. Intrinsic Functions ..
  110:       INTRINSIC          MAX, SQRT
  111: *     ..
  112: *     .. Executable Statements ..
  113: *
  114: *     Test the input parameters.
  115: *
  116:       INFO = 0
  117: *
  118:       IF( LSAME( COMPZ, 'N' ) ) THEN
  119:          ICOMPZ = 0
  120:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
  121:          ICOMPZ = 1
  122:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
  123:          ICOMPZ = 2
  124:       ELSE
  125:          ICOMPZ = -1
  126:       END IF
  127:       IF( ICOMPZ.LT.0 ) THEN
  128:          INFO = -1
  129:       ELSE IF( N.LT.0 ) THEN
  130:          INFO = -2
  131:       ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
  132:      $         N ) ) ) THEN
  133:          INFO = -6
  134:       END IF
  135:       IF( INFO.NE.0 ) THEN
  136:          CALL XERBLA( 'DPTEQR', -INFO )
  137:          RETURN
  138:       END IF
  139: *
  140: *     Quick return if possible
  141: *
  142:       IF( N.EQ.0 )
  143:      $   RETURN
  144: *
  145:       IF( N.EQ.1 ) THEN
  146:          IF( ICOMPZ.GT.0 )
  147:      $      Z( 1, 1 ) = ONE
  148:          RETURN
  149:       END IF
  150:       IF( ICOMPZ.EQ.2 )
  151:      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
  152: *
  153: *     Call DPTTRF to factor the matrix.
  154: *
  155:       CALL DPTTRF( N, D, E, INFO )
  156:       IF( INFO.NE.0 )
  157:      $   RETURN
  158:       DO 10 I = 1, N
  159:          D( I ) = SQRT( D( I ) )
  160:    10 CONTINUE
  161:       DO 20 I = 1, N - 1
  162:          E( I ) = E( I )*D( I )
  163:    20 CONTINUE
  164: *
  165: *     Call DBDSQR to compute the singular values/vectors of the
  166: *     bidiagonal factor.
  167: *
  168:       IF( ICOMPZ.GT.0 ) THEN
  169:          NRU = N
  170:       ELSE
  171:          NRU = 0
  172:       END IF
  173:       CALL DBDSQR( 'Lower', N, 0, NRU, 0, D, E, VT, 1, Z, LDZ, C, 1,
  174:      $             WORK, INFO )
  175: *
  176: *     Square the singular values.
  177: *
  178:       IF( INFO.EQ.0 ) THEN
  179:          DO 30 I = 1, N
  180:             D( I ) = D( I )*D( I )
  181:    30    CONTINUE
  182:       ELSE
  183:          INFO = N + INFO
  184:       END IF
  185: *
  186:       RETURN
  187: *
  188: *     End of DPTEQR
  189: *
  190:       END