File: RELATIVE:/../../../mfix.git/model/DGTSV.f

1     !-----------------------------------------------
2           SUBROUTINE DGTSV( N, NRHS, DL, D, DU, B, LDB, INFO )
3     !-----------------------------------------------
4     !
5     !  -- LAPACK routine (version 3.0) --
6     !     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
7     !     Courant Institute, Argonne National Lab, and Rice University
8     !     October 31, 1999
9     !
10     !     .. Scalar Arguments ..
11           INTEGER            INFO, LDB, N, NRHS
12     !     ..
13     !     .. Array Arguments ..
14           DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * )
15     !     ..
16     !
17     !  Purpose
18     !  =======
19     !
20     !  DGTSV  solves the equation
21     !
22     !     A*X = B,
23     !
24     !  where A is an n by n tridiagonal matrix, by Gaussian elimination with
25     !  partial pivoting.
26     !
27     !  Note that the equation  A'*X = B  may be solved by interchanging the
28     !  order of the arguments DU and DL.
29     !
30     !  Arguments
31     !  =========
32     !
33     !  N       (input) INTEGER
34     !          The order of the matrix A.  N >= 0.
35     !
36     !  NRHS    (input) INTEGER
37     !          The number of right hand sides, i.e., the number of columns
38     !          of the matrix B.  NRHS >= 0.
39     !
40     !  DL      (input/output) DOUBLE PRECISION array, dimension (N-1)
41     !          On entry, DL must contain the (n-1) sub-diagonal elements of
42     !          A.
43     !
44     !          On exit, DL is overwritten by the (n-2) elements of the
45     !          second super-diagonal of the upper triangular matrix U from
46     !          the LU factorization of A, in DL(1), ..., DL(n-2).
47     !
48     !  D       (input/output) DOUBLE PRECISION array, dimension (N)
49     !          On entry, D must contain the diagonal elements of A.
50     !
51     !          On exit, D is overwritten by the n diagonal elements of U.
52     !
53     !  DU      (input/output) DOUBLE PRECISION array, dimension (N-1)
54     !          On entry, DU must contain the (n-1) super-diagonal elements
55     !          of A.
56     !
57     !          On exit, DU is overwritten by the (n-1) elements of the first
58     !          super-diagonal of U.
59     !
60     !  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
61     !          On entry, the N by NRHS matrix of right hand side matrix B.
62     !          On exit, if INFO = 0, the N by NRHS solution matrix X.
63     !
64     !  LDB     (input) INTEGER
65     !          The leading dimension of the array B.  LDB >= max(1,N).
66     !
67     !  INFO    (output) INTEGER
68     !          = 0: successful exit
69     !          < 0: if INFO = -i, the i-th argument had an illegal value
70     !          > 0: if INFO = i, U(i,i) is exactly zero, and the solution
71     !               has not been computed.  The factorization has not been
72     !               completed unless i = N.
73     !
74     !  =====================================================================
75     !
76     !     .. Parameters ..
77           DOUBLE PRECISION   ZERO
78           PARAMETER          ( ZERO = 0.0D+0 )
79     !     ..
80     !     .. Local Scalars ..
81           INTEGER            I, J
82           DOUBLE PRECISION   FACT, TEMP
83     !     ..
84     !     .. Intrinsic Functions ..
85           INTRINSIC          ABS, MAX
86     !     ..
87     !     .. External Subroutines ..
88           EXTERNAL           XERBLA
89     !     ..
90     !     .. Executable Statements ..
91     !
92           INFO = 0
93           IF( N.LT.0 ) THEN
94              INFO = -1
95           ELSE IF( NRHS.LT.0 ) THEN
96              INFO = -2
97           ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
98              INFO = -7
99           END IF
100           IF( INFO.NE.0 ) THEN
101              CALL XERBLA( 'DGTSV ', -INFO )
102              RETURN
103           END IF
104     !
105           IF( N.EQ.0 ) RETURN
106     !
107           IF( NRHS.EQ.1 ) THEN
108              DO 10 I = 1, N - 2
109                 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
110     !
111     !              No row interchange required
112     !
113                    IF( D( I ).NE.ZERO ) THEN
114                       FACT = DL( I ) / D( I )
115                       D( I+1 ) = D( I+1 ) - FACT*DU( I )
116                       B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 )
117                    ELSE
118                       INFO = I
119                       RETURN
120                    END IF
121                    DL( I ) = ZERO
122                 ELSE
123     !
124     !              Interchange rows I and I+1
125     !
126                    FACT = D( I ) / DL( I )
127                    D( I ) = DL( I )
128                    TEMP = D( I+1 )
129                    D( I+1 ) = DU( I ) - FACT*TEMP
130                    DL( I ) = DU( I+1 )
131                    DU( I+1 ) = -FACT*DL( I )
132                    DU( I ) = TEMP
133                    TEMP = B( I, 1 )
134                    B( I, 1 ) = B( I+1, 1 )
135                    B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 )
136                 END IF
137        10    CONTINUE
138              IF( N.GT.1 ) THEN
139                 I = N - 1
140                 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
141                    IF( D( I ).NE.ZERO ) THEN
142                       FACT = DL( I ) / D( I )
143                       D( I+1 ) = D( I+1 ) - FACT*DU( I )
144                       B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 )
145                    ELSE
146                       INFO = I
147                       RETURN
148                    END IF
149                 ELSE
150                    FACT = D( I ) / DL( I )
151                    D( I ) = DL( I )
152                    TEMP = D( I+1 )
153                    D( I+1 ) = DU( I ) - FACT*TEMP
154                    DU( I ) = TEMP
155                    TEMP = B( I, 1 )
156                    B( I, 1 ) = B( I+1, 1 )
157                    B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 )
158                 END IF
159              END IF
160              IF( D( N ).EQ.ZERO ) THEN
161                 INFO = N
162                 RETURN
163              END IF
164           ELSE
165              DO 40 I = 1, N - 2
166                 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
167     !
168     !              No row interchange required
169     !
170                    IF( D( I ).NE.ZERO ) THEN
171                       FACT = DL( I ) / D( I )
172                       D( I+1 ) = D( I+1 ) - FACT*DU( I )
173                       DO 20 J = 1, NRHS
174                          B( I+1, J ) = B( I+1, J ) - FACT*B( I, J )
175        20             CONTINUE
176                    ELSE
177                       INFO = I
178                       RETURN
179                    END IF
180                    DL( I ) = ZERO
181                 ELSE
182     !
183     !              Interchange rows I and I+1
184     !
185                    FACT = D( I ) / DL( I )
186                    D( I ) = DL( I )
187                    TEMP = D( I+1 )
188                    D( I+1 ) = DU( I ) - FACT*TEMP
189                    DL( I ) = DU( I+1 )
190                    DU( I+1 ) = -FACT*DL( I )
191                    DU( I ) = TEMP
192                    DO 30 J = 1, NRHS
193                       TEMP = B( I, J )
194                       B( I, J ) = B( I+1, J )
195                       B( I+1, J ) = TEMP - FACT*B( I+1, J )
196        30          CONTINUE
197                 END IF
198        40    CONTINUE
199              IF( N.GT.1 ) THEN
200                 I = N - 1
201                 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
202                    IF( D( I ).NE.ZERO ) THEN
203                       FACT = DL( I ) / D( I )
204                       D( I+1 ) = D( I+1 ) - FACT*DU( I )
205                       DO 50 J = 1, NRHS
206                          B( I+1, J ) = B( I+1, J ) - FACT*B( I, J )
207        50             CONTINUE
208                    ELSE
209                       INFO = I
210                       RETURN
211                    END IF
212                 ELSE
213                    FACT = D( I ) / DL( I )
214                    D( I ) = DL( I )
215                    TEMP = D( I+1 )
216                    D( I+1 ) = DU( I ) - FACT*TEMP
217                    DU( I ) = TEMP
218                    DO 60 J = 1, NRHS
219                       TEMP = B( I, J )
220                       B( I, J ) = B( I+1, J )
221                       B( I+1, J ) = TEMP - FACT*B( I+1, J )
222        60          CONTINUE
223                 END IF
224              END IF
225              IF( D( N ).EQ.ZERO ) THEN
226                 INFO = N
227                 RETURN
228              END IF
229           END IF
230     !
231     !     Back solve with the matrix U from the factorization.
232     !
233           IF( NRHS.LE.2 ) THEN
234              J = 1
235        70    CONTINUE
236              B( N, J ) = B( N, J ) / D( N )
237              IF( N.GT.1 ) &
238                 B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D( N-1 )
239              DO 80 I = N - 2, 1, -1
240                 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )* &
241                             B( I+2, J ) ) / D( I )
242        80    CONTINUE
243              IF( J.LT.NRHS ) THEN
244                 J = J + 1
245                 GO TO 70
246              END IF
247           ELSE
248              DO 100 J = 1, NRHS
249                 B( N, J ) = B( N, J ) / D( N )
250                 IF( N.GT.1 ) &
251                    B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / &
252                                  D( N-1 )
253                 DO 90 I = N - 2, 1, -1
254                    B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )*&
255                                B( I+2, J ) ) / D( I )
256        90       CONTINUE
257       100    CONTINUE
258           END IF
259     !
260           RETURN
261     !
262     !     End of DGTSV
263     !
264           END SUBROUTINE DGTSV
265     
266