SRC/dgtts2.f

   1 *> \brief \b DGTTS2 solves a system of linear equations with a tridiagonal matrix using the LU factorization computed by sgttrf.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGTTS2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgtts2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgtts2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgtts2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            ITRANS, LDB, N, NRHS
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       INTEGER            IPIV( * )
  28 *       DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DGTTS2 solves one of the systems of equations
  38 *>    A*X = B  or  A**T*X = B,
  39 *> with a tridiagonal matrix A using the LU factorization computed
  40 *> by DGTTRF.
  41 *> \endverbatim
  42 *
  43 *  Arguments:
  44 *  ==========
  45 *
  46 *> \param[in] ITRANS
  47 *> \verbatim
  48 *>          ITRANS is INTEGER
  49 *>          Specifies the form of the system of equations.
  50 *>          = 0:  A * X = B  (No transpose)
  51 *>          = 1:  A**T* X = B  (Transpose)
  52 *>          = 2:  A**T* X = B  (Conjugate transpose = Transpose)
  53 *> \endverbatim
  54 *>
  55 *> \param[in] N
  56 *> \verbatim
  57 *>          N is INTEGER
  58 *>          The order of the matrix A.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] NRHS
  62 *> \verbatim
  63 *>          NRHS is INTEGER
  64 *>          The number of right hand sides, i.e., the number of columns
  65 *>          of the matrix B.  NRHS >= 0.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] DL
  69 *> \verbatim
  70 *>          DL is DOUBLE PRECISION array, dimension (N-1)
  71 *>          The (n-1) multipliers that define the matrix L from the
  72 *>          LU factorization of A.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] D
  76 *> \verbatim
  77 *>          D is DOUBLE PRECISION array, dimension (N)
  78 *>          The n diagonal elements of the upper triangular matrix U from
  79 *>          the LU factorization of A.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] DU
  83 *> \verbatim
  84 *>          DU is DOUBLE PRECISION array, dimension (N-1)
  85 *>          The (n-1) elements of the first super-diagonal of U.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] DU2
  89 *> \verbatim
  90 *>          DU2 is DOUBLE PRECISION array, dimension (N-2)
  91 *>          The (n-2) elements of the second super-diagonal of U.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] IPIV
  95 *> \verbatim
  96 *>          IPIV is INTEGER array, dimension (N)
  97 *>          The pivot indices; for 1 <= i <= n, row i of the matrix was
  98 *>          interchanged with row IPIV(i).  IPIV(i) will always be either
  99 *>          i or i+1; IPIV(i) = i indicates a row interchange was not
 100 *>          required.
 101 *> \endverbatim
 102 *>
 103 *> \param[in,out] B
 104 *> \verbatim
 105 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
 106 *>          On entry, the matrix of right hand side vectors B.
 107 *>          On exit, B is overwritten by the solution vectors X.
 108 *> \endverbatim
 109 *>
 110 *> \param[in] LDB
 111 *> \verbatim
 112 *>          LDB is INTEGER
 113 *>          The leading dimension of the array B.  LDB >= max(1,N).
 114 *> \endverbatim
 115 *
 116 *  Authors:
 117 *  ========
 118 *
 119 *> \author Univ. of Tennessee
 120 *> \author Univ. of California Berkeley
 121 *> \author Univ. of Colorado Denver
 122 *> \author NAG Ltd.
 123 *
 124 *> \date September 2012
 125 *
 126 *> \ingroup doubleGTcomputational
 127 *
 128 *  =====================================================================
 129       SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB )
 130 *
 131 *  -- LAPACK computational routine (version 3.4.2) --
 132 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 133 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 134 *     September 2012
 135 *
 136 *     .. Scalar Arguments ..
 137       INTEGER            ITRANS, LDB, N, NRHS
 138 *     ..
 139 *     .. Array Arguments ..
 140       INTEGER            IPIV( * )
 141       DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
 142 *     ..
 143 *
 144 *  =====================================================================
 145 *
 146 *     .. Local Scalars ..
 147       INTEGER            I, IP, J
 148       DOUBLE PRECISION   TEMP
 149 *     ..
 150 *     .. Executable Statements ..
 151 *
 152 *     Quick return if possible
 153 *
 154       IF( N.EQ.0 .OR. NRHS.EQ.0 )
 155      $   RETURN
 156 *
 157       IF( ITRANS.EQ.0 ) THEN
 158 *
 159 *        Solve A*X = B using the LU factorization of A,
 160 *        overwriting each right hand side vector with its solution.
 161 *
 162          IF( NRHS.LE.1 ) THEN
 163             J = 1
 164    10       CONTINUE
 165 *
 166 *           Solve L*x = b.
 167 *
 168             DO 20 I = 1, N - 1
 169                IP = IPIV( I )
 170                TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J )
 171                B( I, J ) = B( IP, J )
 172                B( I+1, J ) = TEMP
 173    20       CONTINUE
 174 *
 175 *           Solve U*x = b.
 176 *
 177             B( N, J ) = B( N, J ) / D( N )
 178             IF( N.GT.1 )
 179      $         B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
 180      $                       D( N-1 )
 181             DO 30 I = N - 2, 1, -1
 182                B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
 183      $                     B( I+2, J ) ) / D( I )
 184    30       CONTINUE
 185             IF( J.LT.NRHS ) THEN
 186                J = J + 1
 187                GO TO 10
 188             END IF
 189          ELSE
 190             DO 60 J = 1, NRHS
 191 *
 192 *              Solve L*x = b.
 193 *
 194                DO 40 I = 1, N - 1
 195                   IF( IPIV( I ).EQ.I ) THEN
 196                      B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J )
 197                   ELSE
 198                      TEMP = B( I, J )
 199                      B( I, J ) = B( I+1, J )
 200                      B( I+1, J ) = TEMP - DL( I )*B( I, J )
 201                   END IF
 202    40          CONTINUE
 203 *
 204 *              Solve U*x = b.
 205 *
 206                B( N, J ) = B( N, J ) / D( N )
 207                IF( N.GT.1 )
 208      $            B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
 209      $                          D( N-1 )
 210                DO 50 I = N - 2, 1, -1
 211                   B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
 212      $                        B( I+2, J ) ) / D( I )
 213    50          CONTINUE
 214    60       CONTINUE
 215          END IF
 216       ELSE
 217 *
 218 *        Solve A**T * X = B.
 219 *
 220          IF( NRHS.LE.1 ) THEN
 221 *
 222 *           Solve U**T*x = b.
 223 *
 224             J = 1
 225    70       CONTINUE
 226             B( 1, J ) = B( 1, J ) / D( 1 )
 227             IF( N.GT.1 )
 228      $         B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
 229             DO 80 I = 3, N
 230                B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )*
 231      $                     B( I-2, J ) ) / D( I )
 232    80       CONTINUE
 233 *
 234 *           Solve L**T*x = b.
 235 *
 236             DO 90 I = N - 1, 1, -1
 237                IP = IPIV( I )
 238                TEMP = B( I, J ) - DL( I )*B( I+1, J )
 239                B( I, J ) = B( IP, J )
 240                B( IP, J ) = TEMP
 241    90       CONTINUE
 242             IF( J.LT.NRHS ) THEN
 243                J = J + 1
 244                GO TO 70
 245             END IF
 246 *
 247          ELSE
 248             DO 120 J = 1, NRHS
 249 *
 250 *              Solve U**T*x = b.
 251 *
 252                B( 1, J ) = B( 1, J ) / D( 1 )
 253                IF( N.GT.1 )
 254      $            B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
 255                DO 100 I = 3, N
 256                   B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-
 257      $                        DU2( I-2 )*B( I-2, J ) ) / D( I )
 258   100          CONTINUE
 259                DO 110 I = N - 1, 1, -1
 260                   IF( IPIV( I ).EQ.I ) THEN
 261                      B( I, J ) = B( I, J ) - DL( I )*B( I+1, J )
 262                   ELSE
 263                      TEMP = B( I+1, J )
 264                      B( I+1, J ) = B( I, J ) - DL( I )*TEMP
 265                      B( I, J ) = TEMP
 266                   END IF
 267   110          CONTINUE
 268   120       CONTINUE
 269          END IF
 270       END IF
 271 *
 272 *     End of DGTTS2
 273 *
 274       END