1 SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB )
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER ITRANS, LDB, N, NRHS
10 * ..
11 * .. Array Arguments ..
12 INTEGER IPIV( * )
13 DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DGTTS2 solves one of the systems of equations
20 * A*X = B or A**T*X = B,
21 * with a tridiagonal matrix A using the LU factorization computed
22 * by DGTTRF.
23 *
24 * Arguments
25 * =========
26 *
27 * ITRANS (input) INTEGER
28 * Specifies the form of the system of equations.
29 * = 0: A * X = B (No transpose)
30 * = 1: A**T* X = B (Transpose)
31 * = 2: A**T* X = B (Conjugate transpose = Transpose)
32 *
33 * N (input) INTEGER
34 * The order of the matrix A.
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) DOUBLE PRECISION array, dimension (N-1)
41 * The (n-1) multipliers that define the matrix L from the
42 * LU factorization of A.
43 *
44 * D (input) DOUBLE PRECISION array, dimension (N)
45 * The n diagonal elements of the upper triangular matrix U from
46 * the LU factorization of A.
47 *
48 * DU (input) DOUBLE PRECISION array, dimension (N-1)
49 * The (n-1) elements of the first super-diagonal of U.
50 *
51 * DU2 (input) DOUBLE PRECISION array, dimension (N-2)
52 * The (n-2) elements of the second super-diagonal of U.
53 *
54 * IPIV (input) INTEGER array, dimension (N)
55 * The pivot indices; for 1 <= i <= n, row i of the matrix was
56 * interchanged with row IPIV(i). IPIV(i) will always be either
57 * i or i+1; IPIV(i) = i indicates a row interchange was not
58 * required.
59 *
60 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
61 * On entry, the matrix of right hand side vectors B.
62 * On exit, B is overwritten by the solution vectors X.
63 *
64 * LDB (input) INTEGER
65 * The leading dimension of the array B. LDB >= max(1,N).
66 *
67 * =====================================================================
68 *
69 * .. Local Scalars ..
70 INTEGER I, IP, J
71 DOUBLE PRECISION TEMP
72 * ..
73 * .. Executable Statements ..
74 *
75 * Quick return if possible
76 *
77 IF( N.EQ.0 .OR. NRHS.EQ.0 )
78 $ RETURN
79 *
80 IF( ITRANS.EQ.0 ) THEN
81 *
82 * Solve A*X = B using the LU factorization of A,
83 * overwriting each right hand side vector with its solution.
84 *
85 IF( NRHS.LE.1 ) THEN
86 J = 1
87 10 CONTINUE
88 *
89 * Solve L*x = b.
90 *
91 DO 20 I = 1, N - 1
92 IP = IPIV( I )
93 TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J )
94 B( I, J ) = B( IP, J )
95 B( I+1, J ) = TEMP
96 20 CONTINUE
97 *
98 * Solve U*x = b.
99 *
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**T * X = B.
142 *
143 IF( NRHS.LE.1 ) THEN
144 *
145 * Solve U**T*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**T*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**T*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
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER ITRANS, LDB, N, NRHS
10 * ..
11 * .. Array Arguments ..
12 INTEGER IPIV( * )
13 DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DGTTS2 solves one of the systems of equations
20 * A*X = B or A**T*X = B,
21 * with a tridiagonal matrix A using the LU factorization computed
22 * by DGTTRF.
23 *
24 * Arguments
25 * =========
26 *
27 * ITRANS (input) INTEGER
28 * Specifies the form of the system of equations.
29 * = 0: A * X = B (No transpose)
30 * = 1: A**T* X = B (Transpose)
31 * = 2: A**T* X = B (Conjugate transpose = Transpose)
32 *
33 * N (input) INTEGER
34 * The order of the matrix A.
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) DOUBLE PRECISION array, dimension (N-1)
41 * The (n-1) multipliers that define the matrix L from the
42 * LU factorization of A.
43 *
44 * D (input) DOUBLE PRECISION array, dimension (N)
45 * The n diagonal elements of the upper triangular matrix U from
46 * the LU factorization of A.
47 *
48 * DU (input) DOUBLE PRECISION array, dimension (N-1)
49 * The (n-1) elements of the first super-diagonal of U.
50 *
51 * DU2 (input) DOUBLE PRECISION array, dimension (N-2)
52 * The (n-2) elements of the second super-diagonal of U.
53 *
54 * IPIV (input) INTEGER array, dimension (N)
55 * The pivot indices; for 1 <= i <= n, row i of the matrix was
56 * interchanged with row IPIV(i). IPIV(i) will always be either
57 * i or i+1; IPIV(i) = i indicates a row interchange was not
58 * required.
59 *
60 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
61 * On entry, the matrix of right hand side vectors B.
62 * On exit, B is overwritten by the solution vectors X.
63 *
64 * LDB (input) INTEGER
65 * The leading dimension of the array B. LDB >= max(1,N).
66 *
67 * =====================================================================
68 *
69 * .. Local Scalars ..
70 INTEGER I, IP, J
71 DOUBLE PRECISION TEMP
72 * ..
73 * .. Executable Statements ..
74 *
75 * Quick return if possible
76 *
77 IF( N.EQ.0 .OR. NRHS.EQ.0 )
78 $ RETURN
79 *
80 IF( ITRANS.EQ.0 ) THEN
81 *
82 * Solve A*X = B using the LU factorization of A,
83 * overwriting each right hand side vector with its solution.
84 *
85 IF( NRHS.LE.1 ) THEN
86 J = 1
87 10 CONTINUE
88 *
89 * Solve L*x = b.
90 *
91 DO 20 I = 1, N - 1
92 IP = IPIV( I )
93 TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J )
94 B( I, J ) = B( IP, J )
95 B( I+1, J ) = TEMP
96 20 CONTINUE
97 *
98 * Solve U*x = b.
99 *
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**T * X = B.
142 *
143 IF( NRHS.LE.1 ) THEN
144 *
145 * Solve U**T*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**T*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**T*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