001:       SUBROUTINE DGTSV( N, NRHS, DL, D, DU, B, LDB, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            INFO, LDB, N, NRHS
009: *     ..
010: *     .. Array Arguments ..
011:       DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * )
012: *     ..
013: *
014: *  Purpose
015: *  =======
016: *
017: *  DGTSV  solves the equation
018: *
019: *     A*X = B,
020: *
021: *  where A is an n by n tridiagonal matrix, by Gaussian elimination with
022: *  partial pivoting.
023: *
024: *  Note that the equation  A'*X = B  may be solved by interchanging the
025: *  order of the arguments DU and DL.
026: *
027: *  Arguments
028: *  =========
029: *
030: *  N       (input) INTEGER
031: *          The order of the matrix A.  N >= 0.
032: *
033: *  NRHS    (input) INTEGER
034: *          The number of right hand sides, i.e., the number of columns
035: *          of the matrix B.  NRHS >= 0.
036: *
037: *  DL      (input/output) DOUBLE PRECISION array, dimension (N-1)
038: *          On entry, DL must contain the (n-1) sub-diagonal elements of
039: *          A.
040: *
041: *          On exit, DL is overwritten by the (n-2) elements of the
042: *          second super-diagonal of the upper triangular matrix U from
043: *          the LU factorization of A, in DL(1), ..., DL(n-2).
044: *
045: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
046: *          On entry, D must contain the diagonal elements of A.
047: *
048: *          On exit, D is overwritten by the n diagonal elements of U.
049: *
050: *  DU      (input/output) DOUBLE PRECISION array, dimension (N-1)
051: *          On entry, DU must contain the (n-1) super-diagonal elements
052: *          of A.
053: *
054: *          On exit, DU is overwritten by the (n-1) elements of the first
055: *          super-diagonal of U.
056: *
057: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
058: *          On entry, the N by NRHS matrix of right hand side matrix B.
059: *          On exit, if INFO = 0, the N by NRHS solution matrix X.
060: *
061: *  LDB     (input) INTEGER
062: *          The leading dimension of the array B.  LDB >= max(1,N).
063: *
064: *  INFO    (output) INTEGER
065: *          = 0: successful exit
066: *          < 0: if INFO = -i, the i-th argument had an illegal value
067: *          > 0: if INFO = i, U(i,i) is exactly zero, and the solution
068: *               has not been computed.  The factorization has not been
069: *               completed unless i = N.
070: *
071: *  =====================================================================
072: *
073: *     .. Parameters ..
074:       DOUBLE PRECISION   ZERO
075:       PARAMETER          ( ZERO = 0.0D+0 )
076: *     ..
077: *     .. Local Scalars ..
078:       INTEGER            I, J
079:       DOUBLE PRECISION   FACT, TEMP
080: *     ..
081: *     .. Intrinsic Functions ..
082:       INTRINSIC          ABS, MAX
083: *     ..
084: *     .. External Subroutines ..
085:       EXTERNAL           XERBLA
086: *     ..
087: *     .. Executable Statements ..
088: *
089:       INFO = 0
090:       IF( N.LT.0 ) THEN
091:          INFO = -1
092:       ELSE IF( NRHS.LT.0 ) THEN
093:          INFO = -2
094:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
095:          INFO = -7
096:       END IF
097:       IF( INFO.NE.0 ) THEN
098:          CALL XERBLA( 'DGTSV ', -INFO )
099:          RETURN
100:       END IF
101: *
102:       IF( N.EQ.0 )
103:      $   RETURN
104: *
105:       IF( NRHS.EQ.1 ) THEN
106:          DO 10 I = 1, N - 2
107:             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
108: *
109: *              No row interchange required
110: *
111:                IF( D( I ).NE.ZERO ) THEN
112:                   FACT = DL( I ) / D( I )
113:                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
114:                   B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 )
115:                ELSE
116:                   INFO = I
117:                   RETURN
118:                END IF
119:                DL( I ) = ZERO
120:             ELSE
121: *
122: *              Interchange rows I and I+1
123: *
124:                FACT = D( I ) / DL( I )
125:                D( I ) = DL( I )
126:                TEMP = D( I+1 )
127:                D( I+1 ) = DU( I ) - FACT*TEMP
128:                DL( I ) = DU( I+1 )
129:                DU( I+1 ) = -FACT*DL( I )
130:                DU( I ) = TEMP
131:                TEMP = B( I, 1 )
132:                B( I, 1 ) = B( I+1, 1 )
133:                B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 )
134:             END IF
135:    10    CONTINUE
136:          IF( N.GT.1 ) THEN
137:             I = N - 1
138:             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
139:                IF( D( I ).NE.ZERO ) THEN
140:                   FACT = DL( I ) / D( I )
141:                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
142:                   B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 )
143:                ELSE
144:                   INFO = I
145:                   RETURN
146:                END IF
147:             ELSE
148:                FACT = D( I ) / DL( I )
149:                D( I ) = DL( I )
150:                TEMP = D( I+1 )
151:                D( I+1 ) = DU( I ) - FACT*TEMP
152:                DU( I ) = TEMP
153:                TEMP = B( I, 1 )
154:                B( I, 1 ) = B( I+1, 1 )
155:                B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 )
156:             END IF
157:          END IF
158:          IF( D( N ).EQ.ZERO ) THEN
159:             INFO = N
160:             RETURN
161:          END IF
162:       ELSE
163:          DO 40 I = 1, N - 2
164:             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
165: *
166: *              No row interchange required
167: *
168:                IF( D( I ).NE.ZERO ) THEN
169:                   FACT = DL( I ) / D( I )
170:                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
171:                   DO 20 J = 1, NRHS
172:                      B( I+1, J ) = B( I+1, J ) - FACT*B( I, J )
173:    20             CONTINUE
174:                ELSE
175:                   INFO = I
176:                   RETURN
177:                END IF
178:                DL( I ) = ZERO
179:             ELSE
180: *
181: *              Interchange rows I and I+1
182: *
183:                FACT = D( I ) / DL( I )
184:                D( I ) = DL( I )
185:                TEMP = D( I+1 )
186:                D( I+1 ) = DU( I ) - FACT*TEMP
187:                DL( I ) = DU( I+1 )
188:                DU( I+1 ) = -FACT*DL( I )
189:                DU( I ) = TEMP
190:                DO 30 J = 1, NRHS
191:                   TEMP = B( I, J )
192:                   B( I, J ) = B( I+1, J )
193:                   B( I+1, J ) = TEMP - FACT*B( I+1, J )
194:    30          CONTINUE
195:             END IF
196:    40    CONTINUE
197:          IF( N.GT.1 ) THEN
198:             I = N - 1
199:             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
200:                IF( D( I ).NE.ZERO ) THEN
201:                   FACT = DL( I ) / D( I )
202:                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
203:                   DO 50 J = 1, NRHS
204:                      B( I+1, J ) = B( I+1, J ) - FACT*B( I, J )
205:    50             CONTINUE
206:                ELSE
207:                   INFO = I
208:                   RETURN
209:                END IF
210:             ELSE
211:                FACT = D( I ) / DL( I )
212:                D( I ) = DL( I )
213:                TEMP = D( I+1 )
214:                D( I+1 ) = DU( I ) - FACT*TEMP
215:                DU( I ) = TEMP
216:                DO 60 J = 1, NRHS
217:                   TEMP = B( I, J )
218:                   B( I, J ) = B( I+1, J )
219:                   B( I+1, J ) = TEMP - FACT*B( I+1, J )
220:    60          CONTINUE
221:             END IF
222:          END IF
223:          IF( D( N ).EQ.ZERO ) THEN
224:             INFO = N
225:             RETURN
226:          END IF
227:       END IF
228: *
229: *     Back solve with the matrix U from the factorization.
230: *
231:       IF( NRHS.LE.2 ) THEN
232:          J = 1
233:    70    CONTINUE
234:          B( N, J ) = B( N, J ) / D( N )
235:          IF( N.GT.1 )
236:      $      B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D( N-1 )
237:          DO 80 I = N - 2, 1, -1
238:             B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )*
239:      $                  B( I+2, J ) ) / D( I )
240:    80    CONTINUE
241:          IF( J.LT.NRHS ) THEN
242:             J = J + 1
243:             GO TO 70
244:          END IF
245:       ELSE
246:          DO 100 J = 1, NRHS
247:             B( N, J ) = B( N, J ) / D( N )
248:             IF( N.GT.1 )
249:      $         B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
250:      $                       D( N-1 )
251:             DO 90 I = N - 2, 1, -1
252:                B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )*
253:      $                     B( I+2, J ) ) / D( I )
254:    90       CONTINUE
255:   100    CONTINUE
256:       END IF
257: *
258:       RETURN
259: *
260: *     End of DGTSV
261: *
262:       END
263: