Annotation of rpl/lapack/lapack/dgttrf.f, revision 1.1.1.1

1.1       bertrand    1:       SUBROUTINE DGTTRF( N, DL, D, DU, DU2, IPIV, 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:       INTEGER            INFO, N
                     10: *     ..
                     11: *     .. Array Arguments ..
                     12:       INTEGER            IPIV( * )
                     13:       DOUBLE PRECISION   D( * ), DL( * ), DU( * ), DU2( * )
                     14: *     ..
                     15: *
                     16: *  Purpose
                     17: *  =======
                     18: *
                     19: *  DGTTRF computes an LU factorization of a real tridiagonal matrix A
                     20: *  using elimination with partial pivoting and row interchanges.
                     21: *
                     22: *  The factorization has the form
                     23: *     A = L * U
                     24: *  where L is a product of permutation and unit lower bidiagonal
                     25: *  matrices and U is upper triangular with nonzeros in only the main
                     26: *  diagonal and first two superdiagonals.
                     27: *
                     28: *  Arguments
                     29: *  =========
                     30: *
                     31: *  N       (input) INTEGER
                     32: *          The order of the matrix A.
                     33: *
                     34: *  DL      (input/output) DOUBLE PRECISION array, dimension (N-1)
                     35: *          On entry, DL must contain the (n-1) sub-diagonal elements of
                     36: *          A.
                     37: *
                     38: *          On exit, DL is overwritten by the (n-1) multipliers that
                     39: *          define the matrix L from the LU factorization of A.
                     40: *
                     41: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
                     42: *          On entry, D must contain the diagonal elements of A.
                     43: *
                     44: *          On exit, D is overwritten by the n diagonal elements of the
                     45: *          upper triangular matrix U from the LU factorization of A.
                     46: *
                     47: *  DU      (input/output) DOUBLE PRECISION array, dimension (N-1)
                     48: *          On entry, DU must contain the (n-1) super-diagonal elements
                     49: *          of A.
                     50: *
                     51: *          On exit, DU is overwritten by the (n-1) elements of the first
                     52: *          super-diagonal of U.
                     53: *
                     54: *  DU2     (output) DOUBLE PRECISION array, dimension (N-2)
                     55: *          On exit, DU2 is overwritten by the (n-2) elements of the
                     56: *          second super-diagonal of U.
                     57: *
                     58: *  IPIV    (output) INTEGER array, dimension (N)
                     59: *          The pivot indices; for 1 <= i <= n, row i of the matrix was
                     60: *          interchanged with row IPIV(i).  IPIV(i) will always be either
                     61: *          i or i+1; IPIV(i) = i indicates a row interchange was not
                     62: *          required.
                     63: *
                     64: *  INFO    (output) INTEGER
                     65: *          = 0:  successful exit
                     66: *          < 0:  if INFO = -k, the k-th argument had an illegal value
                     67: *          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization
                     68: *                has been completed, but the factor U is exactly
                     69: *                singular, and division by zero will occur if it is used
                     70: *                to solve a system of equations.
                     71: *
                     72: *  =====================================================================
                     73: *
                     74: *     .. Parameters ..
                     75:       DOUBLE PRECISION   ZERO
                     76:       PARAMETER          ( ZERO = 0.0D+0 )
                     77: *     ..
                     78: *     .. Local Scalars ..
                     79:       INTEGER            I
                     80:       DOUBLE PRECISION   FACT, TEMP
                     81: *     ..
                     82: *     .. Intrinsic Functions ..
                     83:       INTRINSIC          ABS
                     84: *     ..
                     85: *     .. External Subroutines ..
                     86:       EXTERNAL           XERBLA
                     87: *     ..
                     88: *     .. Executable Statements ..
                     89: *
                     90:       INFO = 0
                     91:       IF( N.LT.0 ) THEN
                     92:          INFO = -1
                     93:          CALL XERBLA( 'DGTTRF', -INFO )
                     94:          RETURN
                     95:       END IF
                     96: *
                     97: *     Quick return if possible
                     98: *
                     99:       IF( N.EQ.0 )
                    100:      $   RETURN
                    101: *
                    102: *     Initialize IPIV(i) = i and DU2(I) = 0
                    103: *
                    104:       DO 10 I = 1, N
                    105:          IPIV( I ) = I
                    106:    10 CONTINUE
                    107:       DO 20 I = 1, N - 2
                    108:          DU2( I ) = ZERO
                    109:    20 CONTINUE
                    110: *
                    111:       DO 30 I = 1, N - 2
                    112:          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
                    113: *
                    114: *           No row interchange required, eliminate DL(I)
                    115: *
                    116:             IF( D( I ).NE.ZERO ) THEN
                    117:                FACT = DL( I ) / D( I )
                    118:                DL( I ) = FACT
                    119:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
                    120:             END IF
                    121:          ELSE
                    122: *
                    123: *           Interchange rows I and I+1, eliminate DL(I)
                    124: *
                    125:             FACT = D( I ) / DL( I )
                    126:             D( I ) = DL( I )
                    127:             DL( I ) = FACT
                    128:             TEMP = DU( I )
                    129:             DU( I ) = D( I+1 )
                    130:             D( I+1 ) = TEMP - FACT*D( I+1 )
                    131:             DU2( I ) = DU( I+1 )
                    132:             DU( I+1 ) = -FACT*DU( I+1 )
                    133:             IPIV( I ) = I + 1
                    134:          END IF
                    135:    30 CONTINUE
                    136:       IF( N.GT.1 ) THEN
                    137:          I = N - 1
                    138:          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
                    139:             IF( D( I ).NE.ZERO ) THEN
                    140:                FACT = DL( I ) / D( I )
                    141:                DL( I ) = FACT
                    142:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
                    143:             END IF
                    144:          ELSE
                    145:             FACT = D( I ) / DL( I )
                    146:             D( I ) = DL( I )
                    147:             DL( I ) = FACT
                    148:             TEMP = DU( I )
                    149:             DU( I ) = D( I+1 )
                    150:             D( I+1 ) = TEMP - FACT*D( I+1 )
                    151:             IPIV( I ) = I + 1
                    152:          END IF
                    153:       END IF
                    154: *
                    155: *     Check for a zero on the diagonal of U.
                    156: *
                    157:       DO 40 I = 1, N
                    158:          IF( D( I ).EQ.ZERO ) THEN
                    159:             INFO = I
                    160:             GO TO 50
                    161:          END IF
                    162:    40 CONTINUE
                    163:    50 CONTINUE
                    164: *
                    165:       RETURN
                    166: *
                    167: *     End of DGTTRF
                    168: *
                    169:       END

CVSweb interface <joel.bertrand@systella.fr>