Return to dpttrf.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/dpttrf.f, revision 1.8

1.1       bertrand    1:       SUBROUTINE DPTTRF( N, D, E, INFO )
                      2: *
1.8     ! bertrand    3: *  -- LAPACK routine (version 3.3.1) --
1.1       bertrand    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8     ! bertrand    6: *  -- April 2011                                                      --
1.1       bertrand    7: *
                      8: *     .. Scalar Arguments ..
                      9:       INTEGER            INFO, N
                     10: *     ..
                     11: *     .. Array Arguments ..
                     12:       DOUBLE PRECISION   D( * ), E( * )
                     13: *     ..
                     14: *
                     15: *  Purpose
                     16: *  =======
                     17: *
1.8     ! bertrand   18: *  DPTTRF computes the L*D*L**T factorization of a real symmetric
1.1       bertrand   19: *  positive definite tridiagonal matrix A.  The factorization may also
1.8     ! bertrand   20: *  be regarded as having the form A = U**T*D*U.
1.1       bertrand   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
1.8     ! bertrand   31: *          D from the L*D*L**T factorization of A.
1.1       bertrand   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
1.8     ! bertrand   36: *          unit bidiagonal factor L from the L*D*L**T factorization of A.
1.1       bertrand   37: *          E can also be regarded as the superdiagonal of the unit
1.8     ! bertrand   38: *          bidiagonal factor U from the U**T*D*U factorization of A.
1.1       bertrand   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: *
1.8     ! bertrand   80: *     Compute the L*D*L**T (or U**T*D*U) factorization of A.
1.1       bertrand   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