1 /* ./src_f77/zptts2.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
zptts2_(integer * iuplo,integer * n,integer * nrhs,doublereal * d__,doublecomplex * e,doublecomplex * b,integer * ldb)8 /* Subroutine */ int zptts2_(integer *iuplo, integer *n, integer *nrhs,
9 doublereal *d__, doublecomplex *e, doublecomplex *b, integer *ldb)
10 {
11 /* System generated locals */
12 integer b_dim1, b_offset, i__1, i__2, i__3, i__4, i__5, i__6;
13 doublereal d__1;
14 doublecomplex z__1, z__2, z__3, z__4;
15
16 /* Builtin functions */
17 void d_cnjg(doublecomplex *, doublecomplex *);
18
19 /* Local variables */
20 static integer i__, j;
21 extern /* Subroutine */ int zdscal_(integer *, doublereal *,
22 doublecomplex *, integer *);
23
24
25 /* -- LAPACK routine (version 3.0) -- */
26 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
27 /* Courant Institute, Argonne National Lab, and Rice University */
28 /* June 30, 1999 */
29
30 /* .. Scalar Arguments .. */
31 /* .. */
32 /* .. Array Arguments .. */
33 /* .. */
34
35 /* Purpose */
36 /* ======= */
37
38 /* ZPTTS2 solves a tridiagonal system of the form */
39 /* A * X = B */
40 /* using the factorization A = U'*D*U or A = L*D*L' computed by ZPTTRF. */
41 /* D is a diagonal matrix specified in the vector D, U (or L) is a unit */
42 /* bidiagonal matrix whose superdiagonal (subdiagonal) is specified in */
43 /* the vector E, and X and B are N by NRHS matrices. */
44
45 /* Arguments */
46 /* ========= */
47
48 /* IUPLO (input) INTEGER */
49 /* Specifies the form of the factorization and whether the */
50 /* vector E is the superdiagonal of the upper bidiagonal factor */
51 /* U or the subdiagonal of the lower bidiagonal factor L. */
52 /* = 1: A = U'*D*U, E is the superdiagonal of U */
53 /* = 0: A = L*D*L', E is the subdiagonal of L */
54
55 /* N (input) INTEGER */
56 /* The order of the tridiagonal matrix A. N >= 0. */
57
58 /* NRHS (input) INTEGER */
59 /* The number of right hand sides, i.e., the number of columns */
60 /* of the matrix B. NRHS >= 0. */
61
62 /* D (input) DOUBLE PRECISION array, dimension (N) */
63 /* The n diagonal elements of the diagonal matrix D from the */
64 /* factorization A = U'*D*U or A = L*D*L'. */
65
66 /* E (input) COMPLEX*16 array, dimension (N-1) */
67 /* If IUPLO = 1, the (n-1) superdiagonal elements of the unit */
68 /* bidiagonal factor U from the factorization A = U'*D*U. */
69 /* If IUPLO = 0, the (n-1) subdiagonal elements of the unit */
70 /* bidiagonal factor L from the factorization A = L*D*L'. */
71
72 /* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */
73 /* On entry, the right hand side vectors B for the system of */
74 /* linear equations. */
75 /* On exit, the solution vectors, X. */
76
77 /* LDB (input) INTEGER */
78 /* The leading dimension of the array B. LDB >= max(1,N). */
79
80 /* ===================================================================== */
81
82 /* .. Local Scalars .. */
83 /* .. */
84 /* .. External Subroutines .. */
85 /* .. */
86 /* .. Intrinsic Functions .. */
87 /* .. */
88 /* .. Executable Statements .. */
89
90 /* Quick return if possible */
91
92 /* Parameter adjustments */
93 --d__;
94 --e;
95 b_dim1 = *ldb;
96 b_offset = 1 + b_dim1;
97 b -= b_offset;
98
99 /* Function Body */
100 if (*n <= 1) {
101 if (*n == 1) {
102 d__1 = 1. / d__[1];
103 zdscal_(nrhs, &d__1, &b[b_offset], ldb);
104 }
105 return 0;
106 }
107
108 if (*iuplo == 1) {
109
110 /* Solve A * X = B using the factorization A = U'*D*U, */
111 /* overwriting each right hand side vector with its solution. */
112
113 if (*nrhs <= 2) {
114 j = 1;
115 L10:
116
117 /* Solve U' * x = b. */
118
119 i__1 = *n;
120 for (i__ = 2; i__ <= i__1; ++i__) {
121 i__2 = i__ + j * b_dim1;
122 i__3 = i__ + j * b_dim1;
123 i__4 = i__ - 1 + j * b_dim1;
124 d_cnjg(&z__3, &e[i__ - 1]);
125 z__2.r = b[i__4].r * z__3.r - b[i__4].i * z__3.i, z__2.i = b[
126 i__4].r * z__3.i + b[i__4].i * z__3.r;
127 z__1.r = b[i__3].r - z__2.r, z__1.i = b[i__3].i - z__2.i;
128 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
129 /* L20: */
130 }
131
132 /* Solve D * U * x = b. */
133
134 i__1 = *n;
135 for (i__ = 1; i__ <= i__1; ++i__) {
136 i__2 = i__ + j * b_dim1;
137 i__3 = i__ + j * b_dim1;
138 i__4 = i__;
139 z__1.r = b[i__3].r / d__[i__4], z__1.i = b[i__3].i / d__[i__4]
140 ;
141 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
142 /* L30: */
143 }
144 for (i__ = *n - 1; i__ >= 1; --i__) {
145 i__1 = i__ + j * b_dim1;
146 i__2 = i__ + j * b_dim1;
147 i__3 = i__ + 1 + j * b_dim1;
148 i__4 = i__;
149 z__2.r = b[i__3].r * e[i__4].r - b[i__3].i * e[i__4].i,
150 z__2.i = b[i__3].r * e[i__4].i + b[i__3].i * e[i__4]
151 .r;
152 z__1.r = b[i__2].r - z__2.r, z__1.i = b[i__2].i - z__2.i;
153 b[i__1].r = z__1.r, b[i__1].i = z__1.i;
154 /* L40: */
155 }
156 if (j < *nrhs) {
157 ++j;
158 goto L10;
159 }
160 } else {
161 i__1 = *nrhs;
162 for (j = 1; j <= i__1; ++j) {
163
164 /* Solve U' * x = b. */
165
166 i__2 = *n;
167 for (i__ = 2; i__ <= i__2; ++i__) {
168 i__3 = i__ + j * b_dim1;
169 i__4 = i__ + j * b_dim1;
170 i__5 = i__ - 1 + j * b_dim1;
171 d_cnjg(&z__3, &e[i__ - 1]);
172 z__2.r = b[i__5].r * z__3.r - b[i__5].i * z__3.i, z__2.i =
173 b[i__5].r * z__3.i + b[i__5].i * z__3.r;
174 z__1.r = b[i__4].r - z__2.r, z__1.i = b[i__4].i - z__2.i;
175 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
176 /* L50: */
177 }
178
179 /* Solve D * U * x = b. */
180
181 i__2 = *n + j * b_dim1;
182 i__3 = *n + j * b_dim1;
183 i__4 = *n;
184 z__1.r = b[i__3].r / d__[i__4], z__1.i = b[i__3].i / d__[i__4]
185 ;
186 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
187 for (i__ = *n - 1; i__ >= 1; --i__) {
188 i__2 = i__ + j * b_dim1;
189 i__3 = i__ + j * b_dim1;
190 i__4 = i__;
191 z__2.r = b[i__3].r / d__[i__4], z__2.i = b[i__3].i / d__[
192 i__4];
193 i__5 = i__ + 1 + j * b_dim1;
194 i__6 = i__;
195 z__3.r = b[i__5].r * e[i__6].r - b[i__5].i * e[i__6].i,
196 z__3.i = b[i__5].r * e[i__6].i + b[i__5].i * e[
197 i__6].r;
198 z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
199 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
200 /* L60: */
201 }
202 /* L70: */
203 }
204 }
205 } else {
206
207 /* Solve A * X = B using the factorization A = L*D*L', */
208 /* overwriting each right hand side vector with its solution. */
209
210 if (*nrhs <= 2) {
211 j = 1;
212 L80:
213
214 /* Solve L * x = b. */
215
216 i__1 = *n;
217 for (i__ = 2; i__ <= i__1; ++i__) {
218 i__2 = i__ + j * b_dim1;
219 i__3 = i__ + j * b_dim1;
220 i__4 = i__ - 1 + j * b_dim1;
221 i__5 = i__ - 1;
222 z__2.r = b[i__4].r * e[i__5].r - b[i__4].i * e[i__5].i,
223 z__2.i = b[i__4].r * e[i__5].i + b[i__4].i * e[i__5]
224 .r;
225 z__1.r = b[i__3].r - z__2.r, z__1.i = b[i__3].i - z__2.i;
226 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
227 /* L90: */
228 }
229
230 /* Solve D * L' * x = b. */
231
232 i__1 = *n;
233 for (i__ = 1; i__ <= i__1; ++i__) {
234 i__2 = i__ + j * b_dim1;
235 i__3 = i__ + j * b_dim1;
236 i__4 = i__;
237 z__1.r = b[i__3].r / d__[i__4], z__1.i = b[i__3].i / d__[i__4]
238 ;
239 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
240 /* L100: */
241 }
242 for (i__ = *n - 1; i__ >= 1; --i__) {
243 i__1 = i__ + j * b_dim1;
244 i__2 = i__ + j * b_dim1;
245 i__3 = i__ + 1 + j * b_dim1;
246 d_cnjg(&z__3, &e[i__]);
247 z__2.r = b[i__3].r * z__3.r - b[i__3].i * z__3.i, z__2.i = b[
248 i__3].r * z__3.i + b[i__3].i * z__3.r;
249 z__1.r = b[i__2].r - z__2.r, z__1.i = b[i__2].i - z__2.i;
250 b[i__1].r = z__1.r, b[i__1].i = z__1.i;
251 /* L110: */
252 }
253 if (j < *nrhs) {
254 ++j;
255 goto L80;
256 }
257 } else {
258 i__1 = *nrhs;
259 for (j = 1; j <= i__1; ++j) {
260
261 /* Solve L * x = b. */
262
263 i__2 = *n;
264 for (i__ = 2; i__ <= i__2; ++i__) {
265 i__3 = i__ + j * b_dim1;
266 i__4 = i__ + j * b_dim1;
267 i__5 = i__ - 1 + j * b_dim1;
268 i__6 = i__ - 1;
269 z__2.r = b[i__5].r * e[i__6].r - b[i__5].i * e[i__6].i,
270 z__2.i = b[i__5].r * e[i__6].i + b[i__5].i * e[
271 i__6].r;
272 z__1.r = b[i__4].r - z__2.r, z__1.i = b[i__4].i - z__2.i;
273 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
274 /* L120: */
275 }
276
277 /* Solve D * L' * x = b. */
278
279 i__2 = *n + j * b_dim1;
280 i__3 = *n + j * b_dim1;
281 i__4 = *n;
282 z__1.r = b[i__3].r / d__[i__4], z__1.i = b[i__3].i / d__[i__4]
283 ;
284 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
285 for (i__ = *n - 1; i__ >= 1; --i__) {
286 i__2 = i__ + j * b_dim1;
287 i__3 = i__ + j * b_dim1;
288 i__4 = i__;
289 z__2.r = b[i__3].r / d__[i__4], z__2.i = b[i__3].i / d__[
290 i__4];
291 i__5 = i__ + 1 + j * b_dim1;
292 d_cnjg(&z__4, &e[i__]);
293 z__3.r = b[i__5].r * z__4.r - b[i__5].i * z__4.i, z__3.i =
294 b[i__5].r * z__4.i + b[i__5].i * z__4.r;
295 z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
296 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
297 /* L130: */
298 }
299 /* L140: */
300 }
301 }
302 }
303
304 return 0;
305
306 /* End of ZPTTS2 */
307
308 } /* zptts2_ */
309
310