Return to dpttrf.f CVS log

Up to [local] / rpl / lapack / lapack

En route vers la 4.4.1.

    1:       SUBROUTINE DPTTRF( N, D, E, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.3.1) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *  -- April 2011                                                      --
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            INFO, N
   10: *     ..
   11: *     .. Array Arguments ..
   12:       DOUBLE PRECISION   D( * ), E( * )
   13: *     ..
   14: *
   15: *  Purpose
   16: *  =======
   17: *
   18: *  DPTTRF computes the L*D*L**T factorization of a real symmetric
   19: *  positive definite tridiagonal matrix A.  The factorization may also
   20: *  be regarded as having the form A = U**T*D*U.
   21: *
   22: *  Arguments
   23: *  =========
   24: *
   25: *  N       (input) INTEGER
   26: *          The order of the matrix A.  N >= 0.
   27: *
   28: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
   29: *          On entry, the n diagonal elements of the tridiagonal matrix
   30: *          A.  On exit, the n diagonal elements of the diagonal matrix
   31: *          D from the L*D*L**T factorization of A.
   32: *
   33: *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
   34: *          On entry, the (n-1) subdiagonal elements of the tridiagonal
   35: *          matrix A.  On exit, the (n-1) subdiagonal elements of the
   36: *          unit bidiagonal factor L from the L*D*L**T factorization of A.
   37: *          E can also be regarded as the superdiagonal of the unit
   38: *          bidiagonal factor U from the U**T*D*U factorization of A.
   39: *
   40: *  INFO    (output) INTEGER
   41: *          = 0: successful exit
   42: *          < 0: if INFO = -k, the k-th argument had an illegal value
   43: *          > 0: if INFO = k, the leading minor of order k is not
   44: *               positive definite; if k < N, the factorization could not
   45: *               be completed, while if k = N, the factorization was
   46: *               completed, but D(N) <= 0.
   47: *
   48: *  =====================================================================
   49: *
   50: *     .. Parameters ..
   51:       DOUBLE PRECISION   ZERO
   52:       PARAMETER          ( ZERO = 0.0D+0 )
   53: *     ..
   54: *     .. Local Scalars ..
   55:       INTEGER            I, I4
   56:       DOUBLE PRECISION   EI
   57: *     ..
   58: *     .. External Subroutines ..
   59:       EXTERNAL           XERBLA
   60: *     ..
   61: *     .. Intrinsic Functions ..
   62:       INTRINSIC          MOD
   63: *     ..
   64: *     .. Executable Statements ..
   65: *
   66: *     Test the input parameters.
   67: *
   68:       INFO = 0
   69:       IF( N.LT.0 ) THEN
   70:          INFO = -1
   71:          CALL XERBLA( 'DPTTRF', -INFO )
   72:          RETURN
   73:       END IF
   74: *
   75: *     Quick return if possible
   76: *
   77:       IF( N.EQ.0 )
   78:      $   RETURN
   79: *
   80: *     Compute the L*D*L**T (or U**T*D*U) factorization of A.
   81: *
   82:       I4 = MOD( N-1, 4 )
   83:       DO 10 I = 1, I4
   84:          IF( D( I ).LE.ZERO ) THEN
   85:             INFO = I
   86:             GO TO 30
   87:          END IF
   88:          EI = E( I )
   89:          E( I ) = EI / D( I )
   90:          D( I+1 ) = D( I+1 ) - E( I )*EI
   91:    10 CONTINUE
   92: *
   93:       DO 20 I = I4 + 1, N - 4, 4
   94: *
   95: *        Drop out of the loop if d(i) <= 0: the matrix is not positive
   96: *        definite.
   97: *
   98:          IF( D( I ).LE.ZERO ) THEN
   99:             INFO = I
  100:             GO TO 30
  101:          END IF
  102: *
  103: *        Solve for e(i) and d(i+1).
  104: *
  105:          EI = E( I )
  106:          E( I ) = EI / D( I )
  107:          D( I+1 ) = D( I+1 ) - E( I )*EI
  108: *
  109:          IF( D( I+1 ).LE.ZERO ) THEN
  110:             INFO = I + 1
  111:             GO TO 30
  112:          END IF
  113: *
  114: *        Solve for e(i+1) and d(i+2).
  115: *
  116:          EI = E( I+1 )
  117:          E( I+1 ) = EI / D( I+1 )
  118:          D( I+2 ) = D( I+2 ) - E( I+1 )*EI
  119: *
  120:          IF( D( I+2 ).LE.ZERO ) THEN
  121:             INFO = I + 2
  122:             GO TO 30
  123:          END IF
  124: *
  125: *        Solve for e(i+2) and d(i+3).
  126: *
  127:          EI = E( I+2 )
  128:          E( I+2 ) = EI / D( I+2 )
  129:          D( I+3 ) = D( I+3 ) - E( I+2 )*EI
  130: *
  131:          IF( D( I+3 ).LE.ZERO ) THEN
  132:             INFO = I + 3
  133:             GO TO 30
  134:          END IF
  135: *
  136: *        Solve for e(i+3) and d(i+4).
  137: *
  138:          EI = E( I+3 )
  139:          E( I+3 ) = EI / D( I+3 )
  140:          D( I+4 ) = D( I+4 ) - E( I+3 )*EI
  141:    20 CONTINUE
  142: *
  143: *     Check d(n) for positive definiteness.
  144: *
  145:       IF( D( N ).LE.ZERO )
  146:      $   INFO = N
  147: *
  148:    30 CONTINUE
  149:       RETURN
  150: *
  151: *     End of DPTTRF
  152: *
  153:       END