SRC/dgttrf.f

   1 *> \brief \b DGTTRF
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGTTRF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgttrf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgttrf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgttrf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       INTEGER            IPIV( * )
  28 *       DOUBLE PRECISION   D( * ), DL( * ), DU( * ), DU2( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DGTTRF computes an LU factorization of a real tridiagonal matrix A
  38 *> using elimination with partial pivoting and row interchanges.
  39 *>
  40 *> The factorization has the form
  41 *>    A = L * U
  42 *> where L is a product of permutation and unit lower bidiagonal
  43 *> matrices and U is upper triangular with nonzeros in only the main
  44 *> diagonal and first two superdiagonals.
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] N
  51 *> \verbatim
  52 *>          N is INTEGER
  53 *>          The order of the matrix A.
  54 *> \endverbatim
  55 *>
  56 *> \param[in,out] DL
  57 *> \verbatim
  58 *>          DL is DOUBLE PRECISION array, dimension (N-1)
  59 *>          On entry, DL must contain the (n-1) sub-diagonal elements of
  60 *>          A.
  61 *>
  62 *>          On exit, DL is overwritten by the (n-1) multipliers that
  63 *>          define the matrix L from the LU factorization of A.
  64 *> \endverbatim
  65 *>
  66 *> \param[in,out] D
  67 *> \verbatim
  68 *>          D is DOUBLE PRECISION array, dimension (N)
  69 *>          On entry, D must contain the diagonal elements of A.
  70 *>
  71 *>          On exit, D is overwritten by the n diagonal elements of the
  72 *>          upper triangular matrix U from the LU factorization of A.
  73 *> \endverbatim
  74 *>
  75 *> \param[in,out] DU
  76 *> \verbatim
  77 *>          DU is DOUBLE PRECISION array, dimension (N-1)
  78 *>          On entry, DU must contain the (n-1) super-diagonal elements
  79 *>          of A.
  80 *>
  81 *>          On exit, DU is overwritten by the (n-1) elements of the first
  82 *>          super-diagonal of U.
  83 *> \endverbatim
  84 *>
  85 *> \param[out] DU2
  86 *> \verbatim
  87 *>          DU2 is DOUBLE PRECISION array, dimension (N-2)
  88 *>          On exit, DU2 is overwritten by the (n-2) elements of the
  89 *>          second super-diagonal of U.
  90 *> \endverbatim
  91 *>
  92 *> \param[out] IPIV
  93 *> \verbatim
  94 *>          IPIV is INTEGER array, dimension (N)
  95 *>          The pivot indices; for 1 <= i <= n, row i of the matrix was
  96 *>          interchanged with row IPIV(i).  IPIV(i) will always be either
  97 *>          i or i+1; IPIV(i) = i indicates a row interchange was not
  98 *>          required.
  99 *> \endverbatim
 100 *>
 101 *> \param[out] INFO
 102 *> \verbatim
 103 *>          INFO is INTEGER
 104 *>          = 0:  successful exit
 105 *>          < 0:  if INFO = -k, the k-th argument had an illegal value
 106 *>          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization
 107 *>                has been completed, but the factor U is exactly
 108 *>                singular, and division by zero will occur if it is used
 109 *>                to solve a system of equations.
 110 *> \endverbatim
 111 *
 112 *  Authors:
 113 *  ========
 114 *
 115 *> \author Univ. of Tennessee
 116 *> \author Univ. of California Berkeley
 117 *> \author Univ. of Colorado Denver
 118 *> \author NAG Ltd.
 119 *
 120 *> \date September 2012
 121 *
 122 *> \ingroup doubleGTcomputational
 123 *
 124 *  =====================================================================
 125       SUBROUTINE DGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
 126 *
 127 *  -- LAPACK computational routine (version 3.4.2) --
 128 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 129 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 130 *     September 2012
 131 *
 132 *     .. Scalar Arguments ..
 133       INTEGER            INFO, N
 134 *     ..
 135 *     .. Array Arguments ..
 136       INTEGER            IPIV( * )
 137       DOUBLE PRECISION   D( * ), DL( * ), DU( * ), DU2( * )
 138 *     ..
 139 *
 140 *  =====================================================================
 141 *
 142 *     .. Parameters ..
 143       DOUBLE PRECISION   ZERO
 144       PARAMETER          ( ZERO = 0.0D+0 )
 145 *     ..
 146 *     .. Local Scalars ..
 147       INTEGER            I
 148       DOUBLE PRECISION   FACT, TEMP
 149 *     ..
 150 *     .. Intrinsic Functions ..
 151       INTRINSIC          ABS
 152 *     ..
 153 *     .. External Subroutines ..
 154       EXTERNAL           XERBLA
 155 *     ..
 156 *     .. Executable Statements ..
 157 *
 158       INFO = 0
 159       IF( N.LT.0 ) THEN
 160          INFO = -1
 161          CALL XERBLA( 'DGTTRF', -INFO )
 162          RETURN
 163       END IF
 164 *
 165 *     Quick return if possible
 166 *
 167       IF( N.EQ.0 )
 168      $   RETURN
 169 *
 170 *     Initialize IPIV(i) = i and DU2(I) = 0
 171 *
 172       DO 10 I = 1, N
 173          IPIV( I ) = I
 174    10 CONTINUE
 175       DO 20 I = 1, N - 2
 176          DU2( I ) = ZERO
 177    20 CONTINUE
 178 *
 179       DO 30 I = 1, N - 2
 180          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
 181 *
 182 *           No row interchange required, eliminate DL(I)
 183 *
 184             IF( D( I ).NE.ZERO ) THEN
 185                FACT = DL( I ) / D( I )
 186                DL( I ) = FACT
 187                D( I+1 ) = D( I+1 ) - FACT*DU( I )
 188             END IF
 189          ELSE
 190 *
 191 *           Interchange rows I and I+1, eliminate DL(I)
 192 *
 193             FACT = D( I ) / DL( I )
 194             D( I ) = DL( I )
 195             DL( I ) = FACT
 196             TEMP = DU( I )
 197             DU( I ) = D( I+1 )
 198             D( I+1 ) = TEMP - FACT*D( I+1 )
 199             DU2( I ) = DU( I+1 )
 200             DU( I+1 ) = -FACT*DU( I+1 )
 201             IPIV( I ) = I + 1
 202          END IF
 203    30 CONTINUE
 204       IF( N.GT.1 ) THEN
 205          I = N - 1
 206          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
 207             IF( D( I ).NE.ZERO ) THEN
 208                FACT = DL( I ) / D( I )
 209                DL( I ) = FACT
 210                D( I+1 ) = D( I+1 ) - FACT*DU( I )
 211             END IF
 212          ELSE
 213             FACT = D( I ) / DL( I )
 214             D( I ) = DL( I )
 215             DL( I ) = FACT
 216             TEMP = DU( I )
 217             DU( I ) = D( I+1 )
 218             D( I+1 ) = TEMP - FACT*D( I+1 )
 219             IPIV( I ) = I + 1
 220          END IF
 221       END IF
 222 *
 223 *     Check for a zero on the diagonal of U.
 224 *
 225       DO 40 I = 1, N
 226          IF( D( I ).EQ.ZERO ) THEN
 227             INFO = I
 228             GO TO 50
 229          END IF
 230    40 CONTINUE
 231    50 CONTINUE
 232 *
 233       RETURN
 234 *
 235 *     End of DGTTRF
 236 *
 237       END