Mercurial > octave-nkf

comparison libcruft/lapack/zgttrs.f @ 5164:57077d0ddc8e

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