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: *> \ingroup doubleGTcomputational
125: *
126: * =====================================================================
127: SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB )
128: *
129: * -- LAPACK computational routine --
130: * -- LAPACK is a software package provided by Univ. of Tennessee, --
131: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
132: *
133: * .. Scalar Arguments ..
134: INTEGER ITRANS, LDB, N, NRHS
135: * ..
136: * .. Array Arguments ..
137: INTEGER IPIV( * )
138: DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
139: * ..
140: *
141: * =====================================================================
142: *
143: * .. Local Scalars ..
144: INTEGER I, IP, J
145: DOUBLE PRECISION TEMP
146: * ..
147: * .. Executable Statements ..
148: *
149: * Quick return if possible
150: *
151: IF( N.EQ.0 .OR. NRHS.EQ.0 )
152: $ RETURN
153: *
154: IF( ITRANS.EQ.0 ) THEN
155: *
156: * Solve A*X = B using the LU factorization of A,
157: * overwriting each right hand side vector with its solution.
158: *
159: IF( NRHS.LE.1 ) THEN
160: J = 1
161: 10 CONTINUE
162: *
163: * Solve L*x = b.
164: *
165: DO 20 I = 1, N - 1
166: IP = IPIV( I )
167: TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J )
168: B( I, J ) = B( IP, J )
169: B( I+1, J ) = TEMP
170: 20 CONTINUE
171: *
172: * Solve U*x = b.
173: *
174: B( N, J ) = B( N, J ) / D( N )
175: IF( N.GT.1 )
176: $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
177: $ D( N-1 )
178: DO 30 I = N - 2, 1, -1
179: B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
180: $ B( I+2, J ) ) / D( I )
181: 30 CONTINUE
182: IF( J.LT.NRHS ) THEN
183: J = J + 1
184: GO TO 10
185: END IF
186: ELSE
187: DO 60 J = 1, NRHS
188: *
189: * Solve L*x = b.
190: *
191: DO 40 I = 1, N - 1
192: IF( IPIV( I ).EQ.I ) THEN
193: B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J )
194: ELSE
195: TEMP = B( I, J )
196: B( I, J ) = B( I+1, J )
197: B( I+1, J ) = TEMP - DL( I )*B( I, J )
198: END IF
199: 40 CONTINUE
200: *
201: * Solve U*x = b.
202: *
203: B( N, J ) = B( N, J ) / D( N )
204: IF( N.GT.1 )
205: $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
206: $ D( N-1 )
207: DO 50 I = N - 2, 1, -1
208: B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
209: $ B( I+2, J ) ) / D( I )
210: 50 CONTINUE
211: 60 CONTINUE
212: END IF
213: ELSE
214: *
215: * Solve A**T * X = B.
216: *
217: IF( NRHS.LE.1 ) THEN
218: *
219: * Solve U**T*x = b.
220: *
221: J = 1
222: 70 CONTINUE
223: B( 1, J ) = B( 1, J ) / D( 1 )
224: IF( N.GT.1 )
225: $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
226: DO 80 I = 3, N
227: B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )*
228: $ B( I-2, J ) ) / D( I )
229: 80 CONTINUE
230: *
231: * Solve L**T*x = b.
232: *
233: DO 90 I = N - 1, 1, -1
234: IP = IPIV( I )
235: TEMP = B( I, J ) - DL( I )*B( I+1, J )
236: B( I, J ) = B( IP, J )
237: B( IP, J ) = TEMP
238: 90 CONTINUE
239: IF( J.LT.NRHS ) THEN
240: J = J + 1
241: GO TO 70
242: END IF
243: *
244: ELSE
245: DO 120 J = 1, NRHS
246: *
247: * Solve U**T*x = b.
248: *
249: B( 1, J ) = B( 1, J ) / D( 1 )
250: IF( N.GT.1 )
251: $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
252: DO 100 I = 3, N
253: B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-
254: $ DU2( I-2 )*B( I-2, J ) ) / D( I )
255: 100 CONTINUE
256: DO 110 I = N - 1, 1, -1
257: IF( IPIV( I ).EQ.I ) THEN
258: B( I, J ) = B( I, J ) - DL( I )*B( I+1, J )
259: ELSE
260: TEMP = B( I+1, J )
261: B( I+1, J ) = B( I, J ) - DL( I )*TEMP
262: B( I, J ) = TEMP
263: END IF
264: 110 CONTINUE
265: 120 CONTINUE
266: END IF
267: END IF
268: *
269: * End of DGTTS2
270: *
271: END
CVSweb interface <joel.bertrand@systella.fr>