001:       SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            ITRANS, LDB, N, NRHS
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IPIV( * )
013:       DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DGTTS2 solves one of the systems of equations
020: *     A*X = B  or  A'*X = B,
021: *  with a tridiagonal matrix A using the LU factorization computed
022: *  by DGTTRF.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  ITRANS  (input) INTEGER
028: *          Specifies the form of the system of equations.
029: *          = 0:  A * X = B  (No transpose)
030: *          = 1:  A'* X = B  (Transpose)
031: *          = 2:  A'* X = B  (Conjugate transpose = Transpose)
032: *
033: *  N       (input) INTEGER
034: *          The order of the matrix A.
035: *
036: *  NRHS    (input) INTEGER
037: *          The number of right hand sides, i.e., the number of columns
038: *          of the matrix B.  NRHS >= 0.
039: *
040: *  DL      (input) DOUBLE PRECISION array, dimension (N-1)
041: *          The (n-1) multipliers that define the matrix L from the
042: *          LU factorization of A.
043: *
044: *  D       (input) DOUBLE PRECISION array, dimension (N)
045: *          The n diagonal elements of the upper triangular matrix U from
046: *          the LU factorization of A.
047: *
048: *  DU      (input) DOUBLE PRECISION array, dimension (N-1)
049: *          The (n-1) elements of the first super-diagonal of U.
050: *
051: *  DU2     (input) DOUBLE PRECISION array, dimension (N-2)
052: *          The (n-2) elements of the second super-diagonal of U.
053: *
054: *  IPIV    (input) INTEGER array, dimension (N)
055: *          The pivot indices; for 1 <= i <= n, row i of the matrix was
056: *          interchanged with row IPIV(i).  IPIV(i) will always be either
057: *          i or i+1; IPIV(i) = i indicates a row interchange was not
058: *          required.
059: *
060: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
061: *          On entry, the matrix of right hand side vectors B.
062: *          On exit, B is overwritten by the solution vectors X.
063: *
064: *  LDB     (input) INTEGER
065: *          The leading dimension of the array B.  LDB >= max(1,N).
066: *
067: *  =====================================================================
068: *
069: *     .. Local Scalars ..
070:       INTEGER            I, IP, J
071:       DOUBLE PRECISION   TEMP
072: *     ..
073: *     .. Executable Statements ..
074: *
075: *     Quick return if possible
076: *
077:       IF( N.EQ.0 .OR. NRHS.EQ.0 )
078:      $   RETURN
079: *
080:       IF( ITRANS.EQ.0 ) THEN
081: *
082: *        Solve A*X = B using the LU factorization of A,
083: *        overwriting each right hand side vector with its solution.
084: *
085:          IF( NRHS.LE.1 ) THEN
086:             J = 1
087:    10       CONTINUE
088: *
089: *           Solve L*x = b.
090: *
091:             DO 20 I = 1, N - 1
092:                IP = IPIV( I )
093:                TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J )
094:                B( I, J ) = B( IP, J )
095:                B( I+1, J ) = TEMP
096:    20       CONTINUE
097: *
098: *           Solve U*x = b.
099: *
100:             B( N, J ) = B( N, J ) / D( N )
101:             IF( N.GT.1 )
102:      $         B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
103:      $                       D( N-1 )
104:             DO 30 I = N - 2, 1, -1
105:                B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
106:      $                     B( I+2, J ) ) / D( I )
107:    30       CONTINUE
108:             IF( J.LT.NRHS ) THEN
109:                J = J + 1
110:                GO TO 10
111:             END IF
112:          ELSE
113:             DO 60 J = 1, NRHS
114: *
115: *              Solve L*x = b.
116: *
117:                DO 40 I = 1, N - 1
118:                   IF( IPIV( I ).EQ.I ) THEN
119:                      B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J )
120:                   ELSE
121:                      TEMP = B( I, J )
122:                      B( I, J ) = B( I+1, J )
123:                      B( I+1, J ) = TEMP - DL( I )*B( I, J )
124:                   END IF
125:    40          CONTINUE
126: *
127: *              Solve U*x = b.
128: *
129:                B( N, J ) = B( N, J ) / D( N )
130:                IF( N.GT.1 )
131:      $            B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
132:      $                          D( N-1 )
133:                DO 50 I = N - 2, 1, -1
134:                   B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
135:      $                        B( I+2, J ) ) / D( I )
136:    50          CONTINUE
137:    60       CONTINUE
138:          END IF
139:       ELSE
140: *
141: *        Solve A' * X = B.
142: *
143:          IF( NRHS.LE.1 ) THEN
144: *
145: *           Solve U'*x = b.
146: *
147:             J = 1
148:    70       CONTINUE
149:             B( 1, J ) = B( 1, J ) / D( 1 )
150:             IF( N.GT.1 )
151:      $         B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
152:             DO 80 I = 3, N
153:                B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )*
154:      $                     B( I-2, J ) ) / D( I )
155:    80       CONTINUE
156: *
157: *           Solve L'*x = b.
158: *
159:             DO 90 I = N - 1, 1, -1
160:                IP = IPIV( I )
161:                TEMP = B( I, J ) - DL( I )*B( I+1, J )
162:                B( I, J ) = B( IP, J )
163:                B( IP, J ) = TEMP
164:    90       CONTINUE
165:             IF( J.LT.NRHS ) THEN
166:                J = J + 1
167:                GO TO 70
168:             END IF
169: *
170:          ELSE
171:             DO 120 J = 1, NRHS
172: *
173: *              Solve U'*x = b.
174: *
175:                B( 1, J ) = B( 1, J ) / D( 1 )
176:                IF( N.GT.1 )
177:      $            B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
178:                DO 100 I = 3, N
179:                   B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-
180:      $                        DU2( I-2 )*B( I-2, J ) ) / D( I )
181:   100          CONTINUE
182:                DO 110 I = N - 1, 1, -1
183:                   IF( IPIV( I ).EQ.I ) THEN
184:                      B( I, J ) = B( I, J ) - DL( I )*B( I+1, J )
185:                   ELSE
186:                      TEMP = B( I+1, J )
187:                      B( I+1, J ) = B( I, J ) - DL( I )*TEMP
188:                      B( I, J ) = TEMP
189:                   END IF
190:   110          CONTINUE
191:   120       CONTINUE
192:          END IF
193:       END IF
194: *
195: *     End of DGTTS2
196: *
197:       END
198: