LAPACK 3.3.1 Linear Algebra PACKage

# dgttrf.f

Go to the documentation of this file.
```00001       SUBROUTINE DGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            IPIV( * )
00013       DOUBLE PRECISION   D( * ), DL( * ), DU( * ), DU2( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DGTTRF computes an LU factorization of a real tridiagonal matrix A
00020 *  using elimination with partial pivoting and row interchanges.
00021 *
00022 *  The factorization has the form
00023 *     A = L * U
00024 *  where L is a product of permutation and unit lower bidiagonal
00025 *  matrices and U is upper triangular with nonzeros in only the main
00026 *  diagonal and first two superdiagonals.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  N       (input) INTEGER
00032 *          The order of the matrix A.
00033 *
00034 *  DL      (input/output) DOUBLE PRECISION array, dimension (N-1)
00035 *          On entry, DL must contain the (n-1) sub-diagonal elements of
00036 *          A.
00037 *
00038 *          On exit, DL is overwritten by the (n-1) multipliers that
00039 *          define the matrix L from the LU factorization of A.
00040 *
00041 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
00042 *          On entry, D must contain the diagonal elements of A.
00043 *
00044 *          On exit, D is overwritten by the n diagonal elements of the
00045 *          upper triangular matrix U from the LU factorization of A.
00046 *
00047 *  DU      (input/output) DOUBLE PRECISION array, dimension (N-1)
00048 *          On entry, DU must contain the (n-1) super-diagonal elements
00049 *          of A.
00050 *
00051 *          On exit, DU is overwritten by the (n-1) elements of the first
00052 *          super-diagonal of U.
00053 *
00054 *  DU2     (output) DOUBLE PRECISION array, dimension (N-2)
00055 *          On exit, DU2 is overwritten by the (n-2) elements of the
00056 *          second super-diagonal of U.
00057 *
00058 *  IPIV    (output) INTEGER array, dimension (N)
00059 *          The pivot indices; for 1 <= i <= n, row i of the matrix was
00060 *          interchanged with row IPIV(i).  IPIV(i) will always be either
00061 *          i or i+1; IPIV(i) = i indicates a row interchange was not
00062 *          required.
00063 *
00064 *  INFO    (output) INTEGER
00065 *          = 0:  successful exit
00066 *          < 0:  if INFO = -k, the k-th argument had an illegal value
00067 *          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization
00068 *                has been completed, but the factor U is exactly
00069 *                singular, and division by zero will occur if it is used
00070 *                to solve a system of equations.
00071 *
00072 *  =====================================================================
00073 *
00074 *     .. Parameters ..
00075       DOUBLE PRECISION   ZERO
00076       PARAMETER          ( ZERO = 0.0D+0 )
00077 *     ..
00078 *     .. Local Scalars ..
00079       INTEGER            I
00080       DOUBLE PRECISION   FACT, TEMP
00081 *     ..
00082 *     .. Intrinsic Functions ..
00083       INTRINSIC          ABS
00084 *     ..
00085 *     .. External Subroutines ..
00086       EXTERNAL           XERBLA
00087 *     ..
00088 *     .. Executable Statements ..
00089 *
00090       INFO = 0
00091       IF( N.LT.0 ) THEN
00092          INFO = -1
00093          CALL XERBLA( 'DGTTRF', -INFO )
00094          RETURN
00095       END IF
00096 *
00097 *     Quick return if possible
00098 *
00099       IF( N.EQ.0 )
00100      \$   RETURN
00101 *
00102 *     Initialize IPIV(i) = i and DU2(I) = 0
00103 *
00104       DO 10 I = 1, N
00105          IPIV( I ) = I
00106    10 CONTINUE
00107       DO 20 I = 1, N - 2
00108          DU2( I ) = ZERO
00109    20 CONTINUE
00110 *
00111       DO 30 I = 1, N - 2
00112          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
00113 *
00114 *           No row interchange required, eliminate DL(I)
00115 *
00116             IF( D( I ).NE.ZERO ) THEN
00117                FACT = DL( I ) / D( I )
00118                DL( I ) = FACT
00119                D( I+1 ) = D( I+1 ) - FACT*DU( I )
00120             END IF
00121          ELSE
00122 *
00123 *           Interchange rows I and I+1, eliminate DL(I)
00124 *
00125             FACT = D( I ) / DL( I )
00126             D( I ) = DL( I )
00127             DL( I ) = FACT
00128             TEMP = DU( I )
00129             DU( I ) = D( I+1 )
00130             D( I+1 ) = TEMP - FACT*D( I+1 )
00131             DU2( I ) = DU( I+1 )
00132             DU( I+1 ) = -FACT*DU( I+1 )
00133             IPIV( I ) = I + 1
00134          END IF
00135    30 CONTINUE
00136       IF( N.GT.1 ) THEN
00137          I = N - 1
00138          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
00139             IF( D( I ).NE.ZERO ) THEN
00140                FACT = DL( I ) / D( I )
00141                DL( I ) = FACT
00142                D( I+1 ) = D( I+1 ) - FACT*DU( I )
00143             END IF
00144          ELSE
00145             FACT = D( I ) / DL( I )
00146             D( I ) = DL( I )
00147             DL( I ) = FACT
00148             TEMP = DU( I )
00149             DU( I ) = D( I+1 )
00150             D( I+1 ) = TEMP - FACT*D( I+1 )
00151             IPIV( I ) = I + 1
00152          END IF
00153       END IF
00154 *
00155 *     Check for a zero on the diagonal of U.
00156 *
00157       DO 40 I = 1, N
00158          IF( D( I ).EQ.ZERO ) THEN
00159             INFO = I
00160             GO TO 50
00161          END IF
00162    40 CONTINUE
00163    50 CONTINUE
00164 *
00165       RETURN
00166 *
00167 *     End of DGTTRF
00168 *
00169       END
```