Mercurial > octave-nkf

comparison libcruft/lapack/dgtts2.f @ 7053:570a382ce556

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-23 23:17:36 by jwe]
author jwe
date Tue, 23 Oct 2007 23:17:36 +0000
parents
children
comparison
equal deleted inserted replaced
7052:ee70ac66041f 7053:570a382ce556
1 SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB )
2 *
3 * -- LAPACK auxiliary routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER ITRANS, LDB, N, NRHS
9 * ..
10 * .. Array Arguments ..
11 INTEGER IPIV( * )
12 DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DGTTS2 solves one of the systems of equations
19 * A*X = B or A'*X = B,
20 * with a tridiagonal matrix A using the LU factorization computed
21 * by DGTTRF.
22 *
23 * Arguments
24 * =========
25 *
26 * ITRANS (input) INTEGER
27 * Specifies the form of the system of equations.
28 * = 0: A * X = B (No transpose)
29 * = 1: A'* X = B (Transpose)
30 * = 2: A'* X = B (Conjugate transpose = Transpose)
31 *
32 * N (input) INTEGER
33 * The order of the matrix A.
34 *
35 * NRHS (input) INTEGER
36 * The number of right hand sides, i.e., the number of columns
37 * of the matrix B. NRHS >= 0.
38 *
39 * DL (input) DOUBLE PRECISION array, dimension (N-1)
40 * The (n-1) multipliers that define the matrix L from the
41 * LU factorization of A.
42 *
43 * D (input) DOUBLE PRECISION array, dimension (N)
44 * The n diagonal elements of the upper triangular matrix U from
45 * the LU factorization of A.
46 *
47 * DU (input) DOUBLE PRECISION array, dimension (N-1)
48 * The (n-1) elements of the first super-diagonal of U.
49 *
50 * DU2 (input) DOUBLE PRECISION array, dimension (N-2)
51 * The (n-2) elements of the second super-diagonal of U.
52 *
53 * IPIV (input) INTEGER array, dimension (N)
54 * The pivot indices; for 1 <= i <= n, row i of the matrix was
55 * interchanged with row IPIV(i). IPIV(i) will always be either
56 * i or i+1; IPIV(i) = i indicates a row interchange was not
57 * required.
58 *
59 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
60 * On entry, the matrix of right hand side vectors B.
61 * On exit, B is overwritten by the solution vectors X.
62 *
63 * LDB (input) INTEGER
64 * The leading dimension of the array B. LDB >= max(1,N).
65 *
66 * =====================================================================
67 *
68 * .. Local Scalars ..
69 INTEGER I, IP, J
70 DOUBLE PRECISION TEMP
71 * ..
72 * .. Executable Statements ..
73 *
74 * Quick return if possible
75 *
76 IF( N.EQ.0 .OR. NRHS.EQ.0 )
77 $ RETURN
78 *
79 IF( ITRANS.EQ.0 ) THEN
80 *
81 * Solve A*X = B using the LU factorization of A,
82 * overwriting each right hand side vector with its solution.
83 *
84 IF( NRHS.LE.1 ) THEN
85 J = 1
86 10 CONTINUE
87 *
88 * Solve L*x = b.
89 *
90 DO 20 I = 1, N - 1
91 IP = IPIV( I )
92 TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J )
93 B( I, J ) = B( IP, J )
94 B( I+1, J ) = TEMP
95 20 CONTINUE
96 *
97 * Solve U*x = b.
98 *
99 B( N, J ) = B( N, J ) / D( N )
100 IF( N.GT.1 )
101 $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
102 $ D( N-1 )
103 DO 30 I = N - 2, 1, -1
104 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
105 $ B( I+2, J ) ) / D( I )
106 30 CONTINUE
107 IF( J.LT.NRHS ) THEN
108 J = J + 1
109 GO TO 10
110 END IF
111 ELSE
112 DO 60 J = 1, NRHS
113 *
114 * Solve L*x = b.
115 *
116 DO 40 I = 1, N - 1
117 IF( IPIV( I ).EQ.I ) THEN
118 B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J )
119 ELSE
120 TEMP = B( I, J )
121 B( I, J ) = B( I+1, J )
122 B( I+1, J ) = TEMP - DL( I )*B( I, J )
123 END IF
124 40 CONTINUE
125 *
126 * Solve U*x = b.
127 *
128 B( N, J ) = B( N, J ) / D( N )
129 IF( N.GT.1 )
130 $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
131 $ D( N-1 )
132 DO 50 I = N - 2, 1, -1
133 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
134 $ B( I+2, J ) ) / D( I )
135 50 CONTINUE
136 60 CONTINUE
137 END IF
138 ELSE
139 *
140 * Solve A' * X = B.
141 *
142 IF( NRHS.LE.1 ) THEN
143 *
144 * Solve U'*x = b.
145 *
146 J = 1
147 70 CONTINUE
148 B( 1, J ) = B( 1, J ) / D( 1 )
149 IF( N.GT.1 )
150 $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
151 DO 80 I = 3, N
152 B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )*
153 $ B( I-2, J ) ) / D( I )
154 80 CONTINUE
155 *
156 * Solve L'*x = b.
157 *
158 DO 90 I = N - 1, 1, -1
159 IP = IPIV( I )
160 TEMP = B( I, J ) - DL( I )*B( I+1, J )
161 B( I, J ) = B( IP, J )
162 B( IP, J ) = TEMP
163 90 CONTINUE
164 IF( J.LT.NRHS ) THEN
165 J = J + 1
166 GO TO 70
167 END IF
168 *
169 ELSE
170 DO 120 J = 1, NRHS
171 *
172 * Solve U'*x = b.
173 *
174 B( 1, J ) = B( 1, J ) / D( 1 )
175 IF( N.GT.1 )
176 $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
177 DO 100 I = 3, N
178 B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-
179 $ DU2( I-2 )*B( I-2, J ) ) / D( I )
180 100 CONTINUE
181 DO 110 I = N - 1, 1, -1
182 IF( IPIV( I ).EQ.I ) THEN
183 B( I, J ) = B( I, J ) - DL( I )*B( I+1, J )
184 ELSE
185 TEMP = B( I+1, J )
186 B( I+1, J ) = B( I, J ) - DL( I )*TEMP
187 B( I, J ) = TEMP
188 END IF
189 110 CONTINUE
190 120 CONTINUE
191 END IF
192 END IF
193 *
194 * End of DGTTS2
195 *
196 END