Return to zpteqr.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE ZPTEQR( 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( * )
   14:       COMPLEX*16         Z( LDZ, * )
   15: *     ..
   16: *
   17: *  Purpose
   18: *  =======
   19: *
   20: *  ZPTEQR computes all eigenvalues and, optionally, eigenvectors of a
   21: *  symmetric positive definite tridiagonal matrix by first factoring the
   22: *  matrix using DPTTRF and then calling ZBDSQR to compute the singular
   23: *  values of the bidiagonal factor.
   24: *
   25: *  This routine computes the eigenvalues of the positive definite
   26: *  tridiagonal matrix to high relative accuracy.  This means that if the
   27: *  eigenvalues range over many orders of magnitude in size, then the
   28: *  small eigenvalues and corresponding eigenvectors will be computed
   29: *  more accurately than, for example, with the standard QR method.
   30: *
   31: *  The eigenvectors of a full or band positive definite Hermitian matrix
   32: *  can also be found if ZHETRD, ZHPTRD, or ZHBTRD has been used to
   33: *  reduce this matrix to tridiagonal form.  (The reduction to
   34: *  tridiagonal form, however, may preclude the possibility of obtaining
   35: *  high relative accuracy in the small eigenvalues of the original
   36: *  matrix, if these eigenvalues range over many orders of magnitude.)
   37: *
   38: *  Arguments
   39: *  =========
   40: *
   41: *  COMPZ   (input) CHARACTER*1
   42: *          = 'N':  Compute eigenvalues only.
   43: *          = 'V':  Compute eigenvectors of original Hermitian
   44: *                  matrix also.  Array Z contains the unitary matrix
   45: *                  used to reduce the original matrix to tridiagonal
   46: *                  form.
   47: *          = 'I':  Compute eigenvectors of tridiagonal matrix also.
   48: *
   49: *  N       (input) INTEGER
   50: *          The order of the matrix.  N >= 0.
   51: *
   52: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
   53: *          On entry, the n diagonal elements of the tridiagonal 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) COMPLEX*16 array, dimension (LDZ, N)
   63: *          On entry, if COMPZ = 'V', the unitary matrix used in the
   64: *          reduction to tridiagonal form.
   65: *          On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
   66: *          original Hermitian 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:       COMPLEX*16         CZERO, CONE
   94:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
   95:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
   96: *     ..
   97: *     .. External Functions ..
   98:       LOGICAL            LSAME
   99:       EXTERNAL           LSAME
  100: *     ..
  101: *     .. External Subroutines ..
  102:       EXTERNAL           DPTTRF, XERBLA, ZBDSQR, ZLASET
  103: *     ..
  104: *     .. Local Arrays ..
  105:       COMPLEX*16         C( 1, 1 ), VT( 1, 1 )
  106: *     ..
  107: *     .. Local Scalars ..
  108:       INTEGER            I, ICOMPZ, NRU
  109: *     ..
  110: *     .. Intrinsic Functions ..
  111:       INTRINSIC          MAX, SQRT
  112: *     ..
  113: *     .. Executable Statements ..
  114: *
  115: *     Test the input parameters.
  116: *
  117:       INFO = 0
  118: *
  119:       IF( LSAME( COMPZ, 'N' ) ) THEN
  120:          ICOMPZ = 0
  121:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
  122:          ICOMPZ = 1
  123:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
  124:          ICOMPZ = 2
  125:       ELSE
  126:          ICOMPZ = -1
  127:       END IF
  128:       IF( ICOMPZ.LT.0 ) THEN
  129:          INFO = -1
  130:       ELSE IF( N.LT.0 ) THEN
  131:          INFO = -2
  132:       ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
  133:      $         N ) ) ) THEN
  134:          INFO = -6
  135:       END IF
  136:       IF( INFO.NE.0 ) THEN
  137:          CALL XERBLA( 'ZPTEQR', -INFO )
  138:          RETURN
  139:       END IF
  140: *
  141: *     Quick return if possible
  142: *
  143:       IF( N.EQ.0 )
  144:      $   RETURN
  145: *
  146:       IF( N.EQ.1 ) THEN
  147:          IF( ICOMPZ.GT.0 )
  148:      $      Z( 1, 1 ) = CONE
  149:          RETURN
  150:       END IF
  151:       IF( ICOMPZ.EQ.2 )
  152:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
  153: *
  154: *     Call DPTTRF to factor the matrix.
  155: *
  156:       CALL DPTTRF( N, D, E, INFO )
  157:       IF( INFO.NE.0 )
  158:      $   RETURN
  159:       DO 10 I = 1, N
  160:          D( I ) = SQRT( D( I ) )
  161:    10 CONTINUE
  162:       DO 20 I = 1, N - 1
  163:          E( I ) = E( I )*D( I )
  164:    20 CONTINUE
  165: *
  166: *     Call ZBDSQR to compute the singular values/vectors of the
  167: *     bidiagonal factor.
  168: *
  169:       IF( ICOMPZ.GT.0 ) THEN
  170:          NRU = N
  171:       ELSE
  172:          NRU = 0
  173:       END IF
  174:       CALL ZBDSQR( 'Lower', N, 0, NRU, 0, D, E, VT, 1, Z, LDZ, C, 1,
  175:      $             WORK, INFO )
  176: *
  177: *     Square the singular values.
  178: *
  179:       IF( INFO.EQ.0 ) THEN
  180:          DO 30 I = 1, N
  181:             D( I ) = D( I )*D( I )
  182:    30    CONTINUE
  183:       ELSE
  184:          INFO = N + INFO
  185:       END IF
  186: *
  187:       RETURN
  188: *
  189: *     End of ZPTEQR
  190: *
  191:       END