SRC/dgtsv.f

   1 *> \brief <b> DGTSV computes the solution to system of linear equations A * X = B for GT matrices <b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGTSV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgtsv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgtsv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgtsv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGTSV( N, NRHS, DL, D, DU, B, LDB, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, LDB, N, NRHS
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> DGTSV  solves the equation
  37 *>
  38 *>    A*X = B,
  39 *>
  40 *> where A is an n by n tridiagonal matrix, by Gaussian elimination with
  41 *> partial pivoting.
  42 *>
  43 *> Note that the equation  A**T*X = B  may be solved by interchanging the
  44 *> order of the arguments DU and DL.
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] N
  51 *> \verbatim
  52 *>          N is INTEGER
  53 *>          The order of the matrix A.  N >= 0.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] NRHS
  57 *> \verbatim
  58 *>          NRHS is INTEGER
  59 *>          The number of right hand sides, i.e., the number of columns
  60 *>          of the matrix B.  NRHS >= 0.
  61 *> \endverbatim
  62 *>
  63 *> \param[in,out] DL
  64 *> \verbatim
  65 *>          DL is DOUBLE PRECISION array, dimension (N-1)
  66 *>          On entry, DL must contain the (n-1) sub-diagonal elements of
  67 *>          A.
  68 *>
  69 *>          On exit, DL is overwritten by the (n-2) elements of the
  70 *>          second super-diagonal of the upper triangular matrix U from
  71 *>          the LU factorization of A, in DL(1), ..., DL(n-2).
  72 *> \endverbatim
  73 *>
  74 *> \param[in,out] D
  75 *> \verbatim
  76 *>          D is DOUBLE PRECISION array, dimension (N)
  77 *>          On entry, D must contain the diagonal elements of A.
  78 *>
  79 *>          On exit, D is overwritten by the n diagonal elements of U.
  80 *> \endverbatim
  81 *>
  82 *> \param[in,out] DU
  83 *> \verbatim
  84 *>          DU is DOUBLE PRECISION array, dimension (N-1)
  85 *>          On entry, DU must contain the (n-1) super-diagonal elements
  86 *>          of A.
  87 *>
  88 *>          On exit, DU is overwritten by the (n-1) elements of the first
  89 *>          super-diagonal of U.
  90 *> \endverbatim
  91 *>
  92 *> \param[in,out] B
  93 *> \verbatim
  94 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
  95 *>          On entry, the N by NRHS matrix of right hand side matrix B.
  96 *>          On exit, if INFO = 0, the N by NRHS solution matrix X.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] LDB
 100 *> \verbatim
 101 *>          LDB is INTEGER
 102 *>          The leading dimension of the array B.  LDB >= max(1,N).
 103 *> \endverbatim
 104 *>
 105 *> \param[out] INFO
 106 *> \verbatim
 107 *>          INFO is INTEGER
 108 *>          = 0: successful exit
 109 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 110 *>          > 0: if INFO = i, U(i,i) is exactly zero, and the solution
 111 *>               has not been computed.  The factorization has not been
 112 *>               completed unless i = N.
 113 *> \endverbatim
 114 *
 115 *  Authors:
 116 *  ========
 117 *
 118 *> \author Univ. of Tennessee
 119 *> \author Univ. of California Berkeley
 120 *> \author Univ. of Colorado Denver
 121 *> \author NAG Ltd.
 122 *
 123 *> \date September 2012
 124 *
 125 *> \ingroup doubleGTsolve
 126 *
 127 *  =====================================================================
 128       SUBROUTINE DGTSV( N, NRHS, DL, D, DU, B, LDB, INFO )
 129 *
 130 *  -- LAPACK driver routine (version 3.4.2) --
 131 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 132 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 133 *     September 2012
 134 *
 135 *     .. Scalar Arguments ..
 136       INTEGER            INFO, LDB, N, NRHS
 137 *     ..
 138 *     .. Array Arguments ..
 139       DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * )
 140 *     ..
 141 *
 142 *  =====================================================================
 143 *
 144 *     .. Parameters ..
 145       DOUBLE PRECISION   ZERO
 146       PARAMETER          ( ZERO = 0.0D+0 )
 147 *     ..
 148 *     .. Local Scalars ..
 149       INTEGER            I, J
 150       DOUBLE PRECISION   FACT, TEMP
 151 *     ..
 152 *     .. Intrinsic Functions ..
 153       INTRINSIC          ABS, MAX
 154 *     ..
 155 *     .. External Subroutines ..
 156       EXTERNAL           XERBLA
 157 *     ..
 158 *     .. Executable Statements ..
 159 *
 160       INFO = 0
 161       IF( N.LT.0 ) THEN
 162          INFO = -1
 163       ELSE IF( NRHS.LT.0 ) THEN
 164          INFO = -2
 165       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 166          INFO = -7
 167       END IF
 168       IF( INFO.NE.0 ) THEN
 169          CALL XERBLA( 'DGTSV ', -INFO )
 170          RETURN
 171       END IF
 172 *
 173       IF( N.EQ.0 )
 174      $   RETURN
 175 *
 176       IF( NRHS.EQ.1 ) THEN
 177          DO 10 I = 1, N - 2
 178             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
 179 *
 180 *              No row interchange required
 181 *
 182                IF( D( I ).NE.ZERO ) THEN
 183                   FACT = DL( I ) / D( I )
 184                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
 185                   B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 )
 186                ELSE
 187                   INFO = I
 188                   RETURN
 189                END IF
 190                DL( I ) = ZERO
 191             ELSE
 192 *
 193 *              Interchange rows I and I+1
 194 *
 195                FACT = D( I ) / DL( I )
 196                D( I ) = DL( I )
 197                TEMP = D( I+1 )
 198                D( I+1 ) = DU( I ) - FACT*TEMP
 199                DL( I ) = DU( I+1 )
 200                DU( I+1 ) = -FACT*DL( I )
 201                DU( I ) = TEMP
 202                TEMP = B( I, 1 )
 203                B( I, 1 ) = B( I+1, 1 )
 204                B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 )
 205             END IF
 206    10    CONTINUE
 207          IF( N.GT.1 ) THEN
 208             I = N - 1
 209             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
 210                IF( D( I ).NE.ZERO ) THEN
 211                   FACT = DL( I ) / D( I )
 212                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
 213                   B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 )
 214                ELSE
 215                   INFO = I
 216                   RETURN
 217                END IF
 218             ELSE
 219                FACT = D( I ) / DL( I )
 220                D( I ) = DL( I )
 221                TEMP = D( I+1 )
 222                D( I+1 ) = DU( I ) - FACT*TEMP
 223                DU( I ) = TEMP
 224                TEMP = B( I, 1 )
 225                B( I, 1 ) = B( I+1, 1 )
 226                B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 )
 227             END IF
 228          END IF
 229          IF( D( N ).EQ.ZERO ) THEN
 230             INFO = N
 231             RETURN
 232          END IF
 233       ELSE
 234          DO 40 I = 1, N - 2
 235             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
 236 *
 237 *              No row interchange required
 238 *
 239                IF( D( I ).NE.ZERO ) THEN
 240                   FACT = DL( I ) / D( I )
 241                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
 242                   DO 20 J = 1, NRHS
 243                      B( I+1, J ) = B( I+1, J ) - FACT*B( I, J )
 244    20             CONTINUE
 245                ELSE
 246                   INFO = I
 247                   RETURN
 248                END IF
 249                DL( I ) = ZERO
 250             ELSE
 251 *
 252 *              Interchange rows I and I+1
 253 *
 254                FACT = D( I ) / DL( I )
 255                D( I ) = DL( I )
 256                TEMP = D( I+1 )
 257                D( I+1 ) = DU( I ) - FACT*TEMP
 258                DL( I ) = DU( I+1 )
 259                DU( I+1 ) = -FACT*DL( I )
 260                DU( I ) = TEMP
 261                DO 30 J = 1, NRHS
 262                   TEMP = B( I, J )
 263                   B( I, J ) = B( I+1, J )
 264                   B( I+1, J ) = TEMP - FACT*B( I+1, J )
 265    30          CONTINUE
 266             END IF
 267    40    CONTINUE
 268          IF( N.GT.1 ) THEN
 269             I = N - 1
 270             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
 271                IF( D( I ).NE.ZERO ) THEN
 272                   FACT = DL( I ) / D( I )
 273                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
 274                   DO 50 J = 1, NRHS
 275                      B( I+1, J ) = B( I+1, J ) - FACT*B( I, J )
 276    50             CONTINUE
 277                ELSE
 278                   INFO = I
 279                   RETURN
 280                END IF
 281             ELSE
 282                FACT = D( I ) / DL( I )
 283                D( I ) = DL( I )
 284                TEMP = D( I+1 )
 285                D( I+1 ) = DU( I ) - FACT*TEMP
 286                DU( I ) = TEMP
 287                DO 60 J = 1, NRHS
 288                   TEMP = B( I, J )
 289                   B( I, J ) = B( I+1, J )
 290                   B( I+1, J ) = TEMP - FACT*B( I+1, J )
 291    60          CONTINUE
 292             END IF
 293          END IF
 294          IF( D( N ).EQ.ZERO ) THEN
 295             INFO = N
 296             RETURN
 297          END IF
 298       END IF
 299 *
 300 *     Back solve with the matrix U from the factorization.
 301 *
 302       IF( NRHS.LE.2 ) THEN
 303          J = 1
 304    70    CONTINUE
 305          B( N, J ) = B( N, J ) / D( N )
 306          IF( N.GT.1 )
 307      $      B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D( N-1 )
 308          DO 80 I = N - 2, 1, -1
 309             B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )*
 310      $                  B( I+2, J ) ) / D( I )
 311    80    CONTINUE
 312          IF( J.LT.NRHS ) THEN
 313             J = J + 1
 314             GO TO 70
 315          END IF
 316       ELSE
 317          DO 100 J = 1, NRHS
 318             B( N, J ) = B( N, J ) / D( N )
 319             IF( N.GT.1 )
 320      $         B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
 321      $                       D( N-1 )
 322             DO 90 I = N - 2, 1, -1
 323                B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )*
 324      $                     B( I+2, J ) ) / D( I )
 325    90       CONTINUE
 326   100    CONTINUE
 327       END IF
 328 *
 329       RETURN
 330 *
 331 *     End of DGTSV
 332 *
 333       END