|
2329
|
1 SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, |
|
|
2 $ CNORM, INFO ) |
|
|
3 * |
|
7034
|
4 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
6 * November 2006 |
|
2329
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 CHARACTER DIAG, NORMIN, TRANS, UPLO |
|
|
10 INTEGER INFO, LDA, N |
|
|
11 DOUBLE PRECISION SCALE |
|
|
12 * .. |
|
|
13 * .. Array Arguments .. |
|
|
14 DOUBLE PRECISION CNORM( * ) |
|
|
15 COMPLEX*16 A( LDA, * ), X( * ) |
|
|
16 * .. |
|
|
17 * |
|
|
18 * Purpose |
|
|
19 * ======= |
|
|
20 * |
|
|
21 * ZLATRS solves one of the triangular systems |
|
|
22 * |
|
|
23 * A * x = s*b, A**T * x = s*b, or A**H * x = s*b, |
|
|
24 * |
|
|
25 * with scaling to prevent overflow. Here A is an upper or lower |
|
|
26 * triangular matrix, A**T denotes the transpose of A, A**H denotes the |
|
|
27 * conjugate transpose of A, x and b are n-element vectors, and s is a |
|
|
28 * scaling factor, usually less than or equal to 1, chosen so that the |
|
|
29 * components of x will be less than the overflow threshold. If the |
|
|
30 * unscaled problem will not cause overflow, the Level 2 BLAS routine |
|
|
31 * ZTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j), |
|
|
32 * then s is set to 0 and a non-trivial solution to A*x = 0 is returned. |
|
|
33 * |
|
|
34 * Arguments |
|
|
35 * ========= |
|
|
36 * |
|
|
37 * UPLO (input) CHARACTER*1 |
|
|
38 * Specifies whether the matrix A is upper or lower triangular. |
|
|
39 * = 'U': Upper triangular |
|
|
40 * = 'L': Lower triangular |
|
|
41 * |
|
|
42 * TRANS (input) CHARACTER*1 |
|
|
43 * Specifies the operation applied to A. |
|
|
44 * = 'N': Solve A * x = s*b (No transpose) |
|
|
45 * = 'T': Solve A**T * x = s*b (Transpose) |
|
|
46 * = 'C': Solve A**H * x = s*b (Conjugate transpose) |
|
|
47 * |
|
|
48 * DIAG (input) CHARACTER*1 |
|
|
49 * Specifies whether or not the matrix A is unit triangular. |
|
|
50 * = 'N': Non-unit triangular |
|
|
51 * = 'U': Unit triangular |
|
|
52 * |
|
|
53 * NORMIN (input) CHARACTER*1 |
|
|
54 * Specifies whether CNORM has been set or not. |
|
|
55 * = 'Y': CNORM contains the column norms on entry |
|
|
56 * = 'N': CNORM is not set on entry. On exit, the norms will |
|
|
57 * be computed and stored in CNORM. |
|
|
58 * |
|
|
59 * N (input) INTEGER |
|
|
60 * The order of the matrix A. N >= 0. |
|
|
61 * |
|
|
62 * A (input) COMPLEX*16 array, dimension (LDA,N) |
|
|
63 * The triangular matrix A. If UPLO = 'U', the leading n by n |
|
|
64 * upper triangular part of the array A contains the upper |
|
|
65 * triangular matrix, and the strictly lower triangular part of |
|
|
66 * A is not referenced. If UPLO = 'L', the leading n by n lower |
|
|
67 * triangular part of the array A contains the lower triangular |
|
|
68 * matrix, and the strictly upper triangular part of A is not |
|
|
69 * referenced. If DIAG = 'U', the diagonal elements of A are |
|
|
70 * also not referenced and are assumed to be 1. |
|
|
71 * |
|
|
72 * LDA (input) INTEGER |
|
|
73 * The leading dimension of the array A. LDA >= max (1,N). |
|
|
74 * |
|
|
75 * X (input/output) COMPLEX*16 array, dimension (N) |
|
|
76 * On entry, the right hand side b of the triangular system. |
|
|
77 * On exit, X is overwritten by the solution vector x. |
|
|
78 * |
|
|
79 * SCALE (output) DOUBLE PRECISION |
|
|
80 * The scaling factor s for the triangular system |
|
|
81 * A * x = s*b, A**T * x = s*b, or A**H * x = s*b. |
|
|
82 * If SCALE = 0, the matrix A is singular or badly scaled, and |
|
|
83 * the vector x is an exact or approximate solution to A*x = 0. |
|
|
84 * |
|
|
85 * CNORM (input or output) DOUBLE PRECISION array, dimension (N) |
|
|
86 * |
|
|
87 * If NORMIN = 'Y', CNORM is an input argument and CNORM(j) |
|
|
88 * contains the norm of the off-diagonal part of the j-th column |
|
|
89 * of A. If TRANS = 'N', CNORM(j) must be greater than or equal |
|
|
90 * to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) |
|
|
91 * must be greater than or equal to the 1-norm. |
|
|
92 * |
|
|
93 * If NORMIN = 'N', CNORM is an output argument and CNORM(j) |
|
|
94 * returns the 1-norm of the offdiagonal part of the j-th column |
|
|
95 * of A. |
|
|
96 * |
|
|
97 * INFO (output) INTEGER |
|
|
98 * = 0: successful exit |
|
|
99 * < 0: if INFO = -k, the k-th argument had an illegal value |
|
|
100 * |
|
|
101 * Further Details |
|
|
102 * ======= ======= |
|
|
103 * |
|
|
104 * A rough bound on x is computed; if that is less than overflow, ZTRSV |
|
|
105 * is called, otherwise, specific code is used which checks for possible |
|
|
106 * overflow or divide-by-zero at every operation. |
|
|
107 * |
|
|
108 * A columnwise scheme is used for solving A*x = b. The basic algorithm |
|
|
109 * if A is lower triangular is |
|
|
110 * |
|
|
111 * x[1:n] := b[1:n] |
|
|
112 * for j = 1, ..., n |
|
|
113 * x(j) := x(j) / A(j,j) |
|
|
114 * x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] |
|
|
115 * end |
|
|
116 * |
|
|
117 * Define bounds on the components of x after j iterations of the loop: |
|
|
118 * M(j) = bound on x[1:j] |
|
|
119 * G(j) = bound on x[j+1:n] |
|
|
120 * Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. |
|
|
121 * |
|
|
122 * Then for iteration j+1 we have |
|
|
123 * M(j+1) <= G(j) / | A(j+1,j+1) | |
|
|
124 * G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | |
|
|
125 * <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) |
|
|
126 * |
|
|
127 * where CNORM(j+1) is greater than or equal to the infinity-norm of |
|
|
128 * column j+1 of A, not counting the diagonal. Hence |
|
|
129 * |
|
|
130 * G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) |
|
|
131 * 1<=i<=j |
|
|
132 * and |
|
|
133 * |
|
|
134 * |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) |
|
|
135 * 1<=i< j |
|
|
136 * |
|
|
137 * Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTRSV if the |
|
|
138 * reciprocal of the largest M(j), j=1,..,n, is larger than |
|
|
139 * max(underflow, 1/overflow). |
|
|
140 * |
|
|
141 * The bound on x(j) is also used to determine when a step in the |
|
|
142 * columnwise method can be performed without fear of overflow. If |
|
|
143 * the computed bound is greater than a large constant, x is scaled to |
|
|
144 * prevent overflow, but if the bound overflows, x is set to 0, x(j) to |
|
|
145 * 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. |
|
|
146 * |
|
|
147 * Similarly, a row-wise scheme is used to solve A**T *x = b or |
|
|
148 * A**H *x = b. The basic algorithm for A upper triangular is |
|
|
149 * |
|
|
150 * for j = 1, ..., n |
|
|
151 * x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) |
|
|
152 * end |
|
|
153 * |
|
|
154 * We simultaneously compute two bounds |
|
|
155 * G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j |
|
|
156 * M(j) = bound on x(i), 1<=i<=j |
|
|
157 * |
|
|
158 * The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we |
|
|
159 * add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. |
|
|
160 * Then the bound on x(j) is |
|
|
161 * |
|
|
162 * M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | |
|
|
163 * |
|
|
164 * <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) |
|
|
165 * 1<=i<=j |
|
|
166 * |
|
|
167 * and we can safely call ZTRSV if 1/M(n) and 1/G(n) are both greater |
|
|
168 * than max(underflow, 1/overflow). |
|
|
169 * |
|
|
170 * ===================================================================== |
|
|
171 * |
|
|
172 * .. Parameters .. |
|
|
173 DOUBLE PRECISION ZERO, HALF, ONE, TWO |
|
|
174 PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0, |
|
|
175 $ TWO = 2.0D+0 ) |
|
|
176 * .. |
|
|
177 * .. Local Scalars .. |
|
|
178 LOGICAL NOTRAN, NOUNIT, UPPER |
|
|
179 INTEGER I, IMAX, J, JFIRST, JINC, JLAST |
|
|
180 DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL, |
|
|
181 $ XBND, XJ, XMAX |
|
|
182 COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM |
|
|
183 * .. |
|
|
184 * .. External Functions .. |
|
|
185 LOGICAL LSAME |
|
|
186 INTEGER IDAMAX, IZAMAX |
|
|
187 DOUBLE PRECISION DLAMCH, DZASUM |
|
|
188 COMPLEX*16 ZDOTC, ZDOTU, ZLADIV |
|
|
189 EXTERNAL LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC, |
|
|
190 $ ZDOTU, ZLADIV |
|
|
191 * .. |
|
|
192 * .. External Subroutines .. |
|
3333
|
193 EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTRSV |
|
2329
|
194 * .. |
|
|
195 * .. Intrinsic Functions .. |
|
|
196 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN |
|
|
197 * .. |
|
|
198 * .. Statement Functions .. |
|
|
199 DOUBLE PRECISION CABS1, CABS2 |
|
|
200 * .. |
|
|
201 * .. Statement Function definitions .. |
|
|
202 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) |
|
|
203 CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) + |
|
|
204 $ ABS( DIMAG( ZDUM ) / 2.D0 ) |
|
|
205 * .. |
|
|
206 * .. Executable Statements .. |
|
|
207 * |
|
|
208 INFO = 0 |
|
|
209 UPPER = LSAME( UPLO, 'U' ) |
|
|
210 NOTRAN = LSAME( TRANS, 'N' ) |
|
|
211 NOUNIT = LSAME( DIAG, 'N' ) |
|
|
212 * |
|
|
213 * Test the input parameters. |
|
|
214 * |
|
|
215 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN |
|
|
216 INFO = -1 |
|
|
217 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. |
|
|
218 $ LSAME( TRANS, 'C' ) ) THEN |
|
|
219 INFO = -2 |
|
|
220 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN |
|
|
221 INFO = -3 |
|
|
222 ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT. |
|
|
223 $ LSAME( NORMIN, 'N' ) ) THEN |
|
|
224 INFO = -4 |
|
|
225 ELSE IF( N.LT.0 ) THEN |
|
|
226 INFO = -5 |
|
|
227 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN |
|
|
228 INFO = -7 |
|
|
229 END IF |
|
|
230 IF( INFO.NE.0 ) THEN |
|
|
231 CALL XERBLA( 'ZLATRS', -INFO ) |
|
|
232 RETURN |
|
|
233 END IF |
|
|
234 * |
|
|
235 * Quick return if possible |
|
|
236 * |
|
|
237 IF( N.EQ.0 ) |
|
|
238 $ RETURN |
|
|
239 * |
|
|
240 * Determine machine dependent parameters to control overflow. |
|
|
241 * |
|
|
242 SMLNUM = DLAMCH( 'Safe minimum' ) |
|
|
243 BIGNUM = ONE / SMLNUM |
|
|
244 CALL DLABAD( SMLNUM, BIGNUM ) |
|
|
245 SMLNUM = SMLNUM / DLAMCH( 'Precision' ) |
|
|
246 BIGNUM = ONE / SMLNUM |
|
|
247 SCALE = ONE |
|
|
248 * |
|
|
249 IF( LSAME( NORMIN, 'N' ) ) THEN |
|
|
250 * |
|
|
251 * Compute the 1-norm of each column, not including the diagonal. |
|
|
252 * |
|
|
253 IF( UPPER ) THEN |
|
|
254 * |
|
|
255 * A is upper triangular. |
|
|
256 * |
|
|
257 DO 10 J = 1, N |
|
|
258 CNORM( J ) = DZASUM( J-1, A( 1, J ), 1 ) |
|
|
259 10 CONTINUE |
|
|
260 ELSE |
|
|
261 * |
|
|
262 * A is lower triangular. |
|
|
263 * |
|
|
264 DO 20 J = 1, N - 1 |
|
|
265 CNORM( J ) = DZASUM( N-J, A( J+1, J ), 1 ) |
|
|
266 20 CONTINUE |
|
|
267 CNORM( N ) = ZERO |
|
|
268 END IF |
|
|
269 END IF |
|
|
270 * |
|
|
271 * Scale the column norms by TSCAL if the maximum element in CNORM is |
|
|
272 * greater than BIGNUM/2. |
|
|
273 * |
|
|
274 IMAX = IDAMAX( N, CNORM, 1 ) |
|
|
275 TMAX = CNORM( IMAX ) |
|
|
276 IF( TMAX.LE.BIGNUM*HALF ) THEN |
|
|
277 TSCAL = ONE |
|
|
278 ELSE |
|
|
279 TSCAL = HALF / ( SMLNUM*TMAX ) |
|
|
280 CALL DSCAL( N, TSCAL, CNORM, 1 ) |
|
|
281 END IF |
|
|
282 * |
|
|
283 * Compute a bound on the computed solution vector to see if the |
|
|
284 * Level 2 BLAS routine ZTRSV can be used. |
|
|
285 * |
|
|
286 XMAX = ZERO |
|
|
287 DO 30 J = 1, N |
|
|
288 XMAX = MAX( XMAX, CABS2( X( J ) ) ) |
|
|
289 30 CONTINUE |
|
|
290 XBND = XMAX |
|
|
291 * |
|
|
292 IF( NOTRAN ) THEN |
|
|
293 * |
|
|
294 * Compute the growth in A * x = b. |
|
|
295 * |
|
|
296 IF( UPPER ) THEN |
|
|
297 JFIRST = N |
|
|
298 JLAST = 1 |
|
|
299 JINC = -1 |
|
|
300 ELSE |
|
|
301 JFIRST = 1 |
|
|
302 JLAST = N |
|
|
303 JINC = 1 |
|
|
304 END IF |
|
|
305 * |
|
|
306 IF( TSCAL.NE.ONE ) THEN |
|
|
307 GROW = ZERO |
|
|
308 GO TO 60 |
|
|
309 END IF |
|
|
310 * |
|
|
311 IF( NOUNIT ) THEN |
|
|
312 * |
|
|
313 * A is non-unit triangular. |
|
|
314 * |
|
|
315 * Compute GROW = 1/G(j) and XBND = 1/M(j). |
|
|
316 * Initially, G(0) = max{x(i), i=1,...,n}. |
|
|
317 * |
|
|
318 GROW = HALF / MAX( XBND, SMLNUM ) |
|
|
319 XBND = GROW |
|
|
320 DO 40 J = JFIRST, JLAST, JINC |
|
|
321 * |
|
|
322 * Exit the loop if the growth factor is too small. |
|
|
323 * |
|
|
324 IF( GROW.LE.SMLNUM ) |
|
|
325 $ GO TO 60 |
|
|
326 * |
|
|
327 TJJS = A( J, J ) |
|
|
328 TJJ = CABS1( TJJS ) |
|
|
329 * |
|
|
330 IF( TJJ.GE.SMLNUM ) THEN |
|
|
331 * |
|
|
332 * M(j) = G(j-1) / abs(A(j,j)) |
|
|
333 * |
|
|
334 XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) |
|
|
335 ELSE |
|
|
336 * |
|
|
337 * M(j) could overflow, set XBND to 0. |
|
|
338 * |
|
|
339 XBND = ZERO |
|
|
340 END IF |
|
|
341 * |
|
|
342 IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN |
|
|
343 * |
|
|
344 * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) |
|
|
345 * |
|
|
346 GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) |
|
|
347 ELSE |
|
|
348 * |
|
|
349 * G(j) could overflow, set GROW to 0. |
|
|
350 * |
|
|
351 GROW = ZERO |
|
|
352 END IF |
|
|
353 40 CONTINUE |
|
|
354 GROW = XBND |
|
|
355 ELSE |
|
|
356 * |
|
|
357 * A is unit triangular. |
|
|
358 * |
|
|
359 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. |
|
|
360 * |
|
|
361 GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) ) |
|
|
362 DO 50 J = JFIRST, JLAST, JINC |
|
|
363 * |
|
|
364 * Exit the loop if the growth factor is too small. |
|
|
365 * |
|
|
366 IF( GROW.LE.SMLNUM ) |
|
|
367 $ GO TO 60 |
|
|
368 * |
|
|
369 * G(j) = G(j-1)*( 1 + CNORM(j) ) |
|
|
370 * |
|
|
371 GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) |
|
|
372 50 CONTINUE |
|
|
373 END IF |
|
|
374 60 CONTINUE |
|
|
375 * |
|
|
376 ELSE |
|
|
377 * |
|
|
378 * Compute the growth in A**T * x = b or A**H * x = b. |
|
|
379 * |
|
|
380 IF( UPPER ) THEN |
|
|
381 JFIRST = 1 |
|
|
382 JLAST = N |
|
|
383 JINC = 1 |
|
|
384 ELSE |
|
|
385 JFIRST = N |
|
|
386 JLAST = 1 |
|
|
387 JINC = -1 |
|
|
388 END IF |
|
|
389 * |
|
|
390 IF( TSCAL.NE.ONE ) THEN |
|
|
391 GROW = ZERO |
|
|
392 GO TO 90 |
|
|
393 END IF |
|
|
394 * |
|
|
395 IF( NOUNIT ) THEN |
|
|
396 * |
|
|
397 * A is non-unit triangular. |
|
|
398 * |
|
|
399 * Compute GROW = 1/G(j) and XBND = 1/M(j). |
|
|
400 * Initially, M(0) = max{x(i), i=1,...,n}. |
|
|
401 * |
|
|
402 GROW = HALF / MAX( XBND, SMLNUM ) |
|
|
403 XBND = GROW |
|
|
404 DO 70 J = JFIRST, JLAST, JINC |
|
|
405 * |
|
|
406 * Exit the loop if the growth factor is too small. |
|
|
407 * |
|
|
408 IF( GROW.LE.SMLNUM ) |
|
|
409 $ GO TO 90 |
|
|
410 * |
|
|
411 * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) |
|
|
412 * |
|
|
413 XJ = ONE + CNORM( J ) |
|
|
414 GROW = MIN( GROW, XBND / XJ ) |
|
|
415 * |
|
|
416 TJJS = A( J, J ) |
|
|
417 TJJ = CABS1( TJJS ) |
|
|
418 * |
|
|
419 IF( TJJ.GE.SMLNUM ) THEN |
|
|
420 * |
|
|
421 * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) |
|
|
422 * |
|
|
423 IF( XJ.GT.TJJ ) |
|
|
424 $ XBND = XBND*( TJJ / XJ ) |
|
|
425 ELSE |
|
|
426 * |
|
|
427 * M(j) could overflow, set XBND to 0. |
|
|
428 * |
|
|
429 XBND = ZERO |
|
|
430 END IF |
|
|
431 70 CONTINUE |
|
|
432 GROW = MIN( GROW, XBND ) |
|
|
433 ELSE |
|
|
434 * |
|
|
435 * A is unit triangular. |
|
|
436 * |
|
|
437 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. |
|
|
438 * |
|
|
439 GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) ) |
|
|
440 DO 80 J = JFIRST, JLAST, JINC |
|
|
441 * |
|
|
442 * Exit the loop if the growth factor is too small. |
|
|
443 * |
|
|
444 IF( GROW.LE.SMLNUM ) |
|
|
445 $ GO TO 90 |
|
|
446 * |
|
|
447 * G(j) = ( 1 + CNORM(j) )*G(j-1) |
|
|
448 * |
|
|
449 XJ = ONE + CNORM( J ) |
|
|
450 GROW = GROW / XJ |
|
|
451 80 CONTINUE |
|
|
452 END IF |
|
|
453 90 CONTINUE |
|
|
454 END IF |
|
|
455 * |
|
|
456 IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN |
|
|
457 * |
|
|
458 * Use the Level 2 BLAS solve if the reciprocal of the bound on |
|
|
459 * elements of X is not too small. |
|
|
460 * |
|
|
461 CALL ZTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 ) |
|
|
462 ELSE |
|
|
463 * |
|
|
464 * Use a Level 1 BLAS solve, scaling intermediate results. |
|
|
465 * |
|
|
466 IF( XMAX.GT.BIGNUM*HALF ) THEN |
|
|
467 * |
|
|
468 * Scale X so that its components are less than or equal to |
|
|
469 * BIGNUM in absolute value. |
|
|
470 * |
|
|
471 SCALE = ( BIGNUM*HALF ) / XMAX |
|
|
472 CALL ZDSCAL( N, SCALE, X, 1 ) |
|
|
473 XMAX = BIGNUM |
|
|
474 ELSE |
|
|
475 XMAX = XMAX*TWO |
|
|
476 END IF |
|
|
477 * |
|
|
478 IF( NOTRAN ) THEN |
|
|
479 * |
|
|
480 * Solve A * x = b |
|
|
481 * |
|
|
482 DO 120 J = JFIRST, JLAST, JINC |
|
|
483 * |
|
|
484 * Compute x(j) = b(j) / A(j,j), scaling x if necessary. |
|
|
485 * |
|
|
486 XJ = CABS1( X( J ) ) |
|
|
487 IF( NOUNIT ) THEN |
|
|
488 TJJS = A( J, J )*TSCAL |
|
|
489 ELSE |
|
|
490 TJJS = TSCAL |
|
|
491 IF( TSCAL.EQ.ONE ) |
|
|
492 $ GO TO 110 |
|
|
493 END IF |
|
|
494 TJJ = CABS1( TJJS ) |
|
|
495 IF( TJJ.GT.SMLNUM ) THEN |
|
|
496 * |
|
|
497 * abs(A(j,j)) > SMLNUM: |
|
|
498 * |
|
|
499 IF( TJJ.LT.ONE ) THEN |
|
|
500 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
501 * |
|
|
502 * Scale x by 1/b(j). |
|
|
503 * |
|
|
504 REC = ONE / XJ |
|
|
505 CALL ZDSCAL( N, REC, X, 1 ) |
|
|
506 SCALE = SCALE*REC |
|
|
507 XMAX = XMAX*REC |
|
|
508 END IF |
|
|
509 END IF |
|
|
510 X( J ) = ZLADIV( X( J ), TJJS ) |
|
|
511 XJ = CABS1( X( J ) ) |
|
|
512 ELSE IF( TJJ.GT.ZERO ) THEN |
|
|
513 * |
|
|
514 * 0 < abs(A(j,j)) <= SMLNUM: |
|
|
515 * |
|
|
516 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
517 * |
|
|
518 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM |
|
|
519 * to avoid overflow when dividing by A(j,j). |
|
|
520 * |
|
|
521 REC = ( TJJ*BIGNUM ) / XJ |
|
|
522 IF( CNORM( J ).GT.ONE ) THEN |
|
|
523 * |
|
|
524 * Scale by 1/CNORM(j) to avoid overflow when |
|
|
525 * multiplying x(j) times column j. |
|
|
526 * |
|
|
527 REC = REC / CNORM( J ) |
|
|
528 END IF |
|
|
529 CALL ZDSCAL( N, REC, X, 1 ) |
|
|
530 SCALE = SCALE*REC |
|
|
531 XMAX = XMAX*REC |
|
|
532 END IF |
|
|
533 X( J ) = ZLADIV( X( J ), TJJS ) |
|
|
534 XJ = CABS1( X( J ) ) |
|
|
535 ELSE |
|
|
536 * |
|
|
537 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and |
|
|
538 * scale = 0, and compute a solution to A*x = 0. |
|
|
539 * |
|
|
540 DO 100 I = 1, N |
|
|
541 X( I ) = ZERO |
|
|
542 100 CONTINUE |
|
|
543 X( J ) = ONE |
|
|
544 XJ = ONE |
|
|
545 SCALE = ZERO |
|
|
546 XMAX = ZERO |
|
|
547 END IF |
|
|
548 110 CONTINUE |
|
|
549 * |
|
|
550 * Scale x if necessary to avoid overflow when adding a |
|
|
551 * multiple of column j of A. |
|
|
552 * |
|
|
553 IF( XJ.GT.ONE ) THEN |
|
|
554 REC = ONE / XJ |
|
|
555 IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN |
|
|
556 * |
|
|
557 * Scale x by 1/(2*abs(x(j))). |
|
|
558 * |
|
|
559 REC = REC*HALF |
|
|
560 CALL ZDSCAL( N, REC, X, 1 ) |
|
|
561 SCALE = SCALE*REC |
|
|
562 END IF |
|
|
563 ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN |
|
|
564 * |
|
|
565 * Scale x by 1/2. |
|
|
566 * |
|
|
567 CALL ZDSCAL( N, HALF, X, 1 ) |
|
|
568 SCALE = SCALE*HALF |
|
|
569 END IF |
|
|
570 * |
|
|
571 IF( UPPER ) THEN |
|
|
572 IF( J.GT.1 ) THEN |
|
|
573 * |
|
|
574 * Compute the update |
|
|
575 * x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) |
|
|
576 * |
|
|
577 CALL ZAXPY( J-1, -X( J )*TSCAL, A( 1, J ), 1, X, |
|
|
578 $ 1 ) |
|
|
579 I = IZAMAX( J-1, X, 1 ) |
|
|
580 XMAX = CABS1( X( I ) ) |
|
|
581 END IF |
|
|
582 ELSE |
|
|
583 IF( J.LT.N ) THEN |
|
|
584 * |
|
|
585 * Compute the update |
|
|
586 * x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) |
|
|
587 * |
|
|
588 CALL ZAXPY( N-J, -X( J )*TSCAL, A( J+1, J ), 1, |
|
|
589 $ X( J+1 ), 1 ) |
|
|
590 I = J + IZAMAX( N-J, X( J+1 ), 1 ) |
|
|
591 XMAX = CABS1( X( I ) ) |
|
|
592 END IF |
|
|
593 END IF |
|
|
594 120 CONTINUE |
|
|
595 * |
|
|
596 ELSE IF( LSAME( TRANS, 'T' ) ) THEN |
|
|
597 * |
|
|
598 * Solve A**T * x = b |
|
|
599 * |
|
|
600 DO 170 J = JFIRST, JLAST, JINC |
|
|
601 * |
|
|
602 * Compute x(j) = b(j) - sum A(k,j)*x(k). |
|
|
603 * k<>j |
|
|
604 * |
|
|
605 XJ = CABS1( X( J ) ) |
|
|
606 USCAL = TSCAL |
|
|
607 REC = ONE / MAX( XMAX, ONE ) |
|
|
608 IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN |
|
|
609 * |
|
|
610 * If x(j) could overflow, scale x by 1/(2*XMAX). |
|
|
611 * |
|
|
612 REC = REC*HALF |
|
|
613 IF( NOUNIT ) THEN |
|
|
614 TJJS = A( J, J )*TSCAL |
|
|
615 ELSE |
|
|
616 TJJS = TSCAL |
|
|
617 END IF |
|
|
618 TJJ = CABS1( TJJS ) |
|
|
619 IF( TJJ.GT.ONE ) THEN |
|
|
620 * |
|
|
621 * Divide by A(j,j) when scaling x if A(j,j) > 1. |
|
|
622 * |
|
|
623 REC = MIN( ONE, REC*TJJ ) |
|
|
624 USCAL = ZLADIV( USCAL, TJJS ) |
|
|
625 END IF |
|
|
626 IF( REC.LT.ONE ) THEN |
|
|
627 CALL ZDSCAL( N, REC, X, 1 ) |
|
|
628 SCALE = SCALE*REC |
|
|
629 XMAX = XMAX*REC |
|
|
630 END IF |
|
|
631 END IF |
|
|
632 * |
|
|
633 CSUMJ = ZERO |
|
|
634 IF( USCAL.EQ.DCMPLX( ONE ) ) THEN |
|
|
635 * |
|
|
636 * If the scaling needed for A in the dot product is 1, |
|
|
637 * call ZDOTU to perform the dot product. |
|
|
638 * |
|
|
639 IF( UPPER ) THEN |
|
|
640 CSUMJ = ZDOTU( J-1, A( 1, J ), 1, X, 1 ) |
|
|
641 ELSE IF( J.LT.N ) THEN |
|
|
642 CSUMJ = ZDOTU( N-J, A( J+1, J ), 1, X( J+1 ), 1 ) |
|
|
643 END IF |
|
|
644 ELSE |
|
|
645 * |
|
|
646 * Otherwise, use in-line code for the dot product. |
|
|
647 * |
|
|
648 IF( UPPER ) THEN |
|
|
649 DO 130 I = 1, J - 1 |
|
|
650 CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I ) |
|
|
651 130 CONTINUE |
|
|
652 ELSE IF( J.LT.N ) THEN |
|
|
653 DO 140 I = J + 1, N |
|
|
654 CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I ) |
|
|
655 140 CONTINUE |
|
|
656 END IF |
|
|
657 END IF |
|
|
658 * |
|
|
659 IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN |
|
|
660 * |
|
|
661 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) |
|
|
662 * was not used to scale the dotproduct. |
|
|
663 * |
|
|
664 X( J ) = X( J ) - CSUMJ |
|
|
665 XJ = CABS1( X( J ) ) |
|
|
666 IF( NOUNIT ) THEN |
|
|
667 TJJS = A( J, J )*TSCAL |
|
|
668 ELSE |
|
|
669 TJJS = TSCAL |
|
|
670 IF( TSCAL.EQ.ONE ) |
|
|
671 $ GO TO 160 |
|
|
672 END IF |
|
|
673 * |
|
|
674 * Compute x(j) = x(j) / A(j,j), scaling if necessary. |
|
|
675 * |
|
|
676 TJJ = CABS1( TJJS ) |
|
|
677 IF( TJJ.GT.SMLNUM ) THEN |
|
|
678 * |
|
|
679 * abs(A(j,j)) > SMLNUM: |
|
|
680 * |
|
|
681 IF( TJJ.LT.ONE ) THEN |
|
|
682 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
683 * |
|
|
684 * Scale X by 1/abs(x(j)). |
|
|
685 * |
|
|
686 REC = ONE / XJ |
|
|
687 CALL ZDSCAL( N, REC, X, 1 ) |
|
|
688 SCALE = SCALE*REC |
|
|
689 XMAX = XMAX*REC |
|
|
690 END IF |
|
|
691 END IF |
|
|
692 X( J ) = ZLADIV( X( J ), TJJS ) |
|
|
693 ELSE IF( TJJ.GT.ZERO ) THEN |
|
|
694 * |
|
|
695 * 0 < abs(A(j,j)) <= SMLNUM: |
|
|
696 * |
|
|
697 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
698 * |
|
|
699 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. |
|
|
700 * |
|
|
701 REC = ( TJJ*BIGNUM ) / XJ |
|
|
702 CALL ZDSCAL( N, REC, X, 1 ) |
|
|
703 SCALE = SCALE*REC |
|
|
704 XMAX = XMAX*REC |
|
|
705 END IF |
|
|
706 X( J ) = ZLADIV( X( J ), TJJS ) |
|
|
707 ELSE |
|
|
708 * |
|
|
709 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and |
|
|
710 * scale = 0 and compute a solution to A**T *x = 0. |
|
|
711 * |
|
|
712 DO 150 I = 1, N |
|
|
713 X( I ) = ZERO |
|
|
714 150 CONTINUE |
|
|
715 X( J ) = ONE |
|
|
716 SCALE = ZERO |
|
|
717 XMAX = ZERO |
|
|
718 END IF |
|
|
719 160 CONTINUE |
|
|
720 ELSE |
|
|
721 * |
|
|
722 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot |
|
|
723 * product has already been divided by 1/A(j,j). |
|
|
724 * |
|
|
725 X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ |
|
|
726 END IF |
|
|
727 XMAX = MAX( XMAX, CABS1( X( J ) ) ) |
|
|
728 170 CONTINUE |
|
|
729 * |
|
|
730 ELSE |
|
|
731 * |
|
|
732 * Solve A**H * x = b |
|
|
733 * |
|
|
734 DO 220 J = JFIRST, JLAST, JINC |
|
|
735 * |
|
|
736 * Compute x(j) = b(j) - sum A(k,j)*x(k). |
|
|
737 * k<>j |
|
|
738 * |
|
|
739 XJ = CABS1( X( J ) ) |
|
|
740 USCAL = TSCAL |
|
|
741 REC = ONE / MAX( XMAX, ONE ) |
|
|
742 IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN |
|
|
743 * |
|
|
744 * If x(j) could overflow, scale x by 1/(2*XMAX). |
|
|
745 * |
|
|
746 REC = REC*HALF |
|
|
747 IF( NOUNIT ) THEN |
|
|
748 TJJS = DCONJG( A( J, J ) )*TSCAL |
|
|
749 ELSE |
|
|
750 TJJS = TSCAL |
|
|
751 END IF |
|
|
752 TJJ = CABS1( TJJS ) |
|
|
753 IF( TJJ.GT.ONE ) THEN |
|
|
754 * |
|
|
755 * Divide by A(j,j) when scaling x if A(j,j) > 1. |
|
|
756 * |
|
|
757 REC = MIN( ONE, REC*TJJ ) |
|
|
758 USCAL = ZLADIV( USCAL, TJJS ) |
|
|
759 END IF |
|
|
760 IF( REC.LT.ONE ) THEN |
|
|
761 CALL ZDSCAL( N, REC, X, 1 ) |
|
|
762 SCALE = SCALE*REC |
|
|
763 XMAX = XMAX*REC |
|
|
764 END IF |
|
|
765 END IF |
|
|
766 * |
|
|
767 CSUMJ = ZERO |
|
|
768 IF( USCAL.EQ.DCMPLX( ONE ) ) THEN |
|
|
769 * |
|
|
770 * If the scaling needed for A in the dot product is 1, |
|
|
771 * call ZDOTC to perform the dot product. |
|
|
772 * |
|
|
773 IF( UPPER ) THEN |
|
|
774 CSUMJ = ZDOTC( J-1, A( 1, J ), 1, X, 1 ) |
|
|
775 ELSE IF( J.LT.N ) THEN |
|
|
776 CSUMJ = ZDOTC( N-J, A( J+1, J ), 1, X( J+1 ), 1 ) |
|
|
777 END IF |
|
|
778 ELSE |
|
|
779 * |
|
|
780 * Otherwise, use in-line code for the dot product. |
|
|
781 * |
|
|
782 IF( UPPER ) THEN |
|
|
783 DO 180 I = 1, J - 1 |
|
|
784 CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )* |
|
|
785 $ X( I ) |
|
|
786 180 CONTINUE |
|
|
787 ELSE IF( J.LT.N ) THEN |
|
|
788 DO 190 I = J + 1, N |
|
|
789 CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )* |
|
|
790 $ X( I ) |
|
|
791 190 CONTINUE |
|
|
792 END IF |
|
|
793 END IF |
|
|
794 * |
|
|
795 IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN |
|
|
796 * |
|
|
797 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) |
|
|
798 * was not used to scale the dotproduct. |
|
|
799 * |
|
|
800 X( J ) = X( J ) - CSUMJ |
|
|
801 XJ = CABS1( X( J ) ) |
|
|
802 IF( NOUNIT ) THEN |
|
|
803 TJJS = DCONJG( A( J, J ) )*TSCAL |
|
|
804 ELSE |
|
|
805 TJJS = TSCAL |
|
|
806 IF( TSCAL.EQ.ONE ) |
|
|
807 $ GO TO 210 |
|
|
808 END IF |
|
|
809 * |
|
|
810 * Compute x(j) = x(j) / A(j,j), scaling if necessary. |
|
|
811 * |
|
|
812 TJJ = CABS1( TJJS ) |
|
|
813 IF( TJJ.GT.SMLNUM ) THEN |
|
|
814 * |
|
|
815 * abs(A(j,j)) > SMLNUM: |
|
|
816 * |
|
|
817 IF( TJJ.LT.ONE ) THEN |
|
|
818 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
819 * |
|
|
820 * Scale X by 1/abs(x(j)). |
|
|
821 * |
|
|
822 REC = ONE / XJ |
|
|
823 CALL ZDSCAL( N, REC, X, 1 ) |
|
|
824 SCALE = SCALE*REC |
|
|
825 XMAX = XMAX*REC |
|
|
826 END IF |
|
|
827 END IF |
|
|
828 X( J ) = ZLADIV( X( J ), TJJS ) |
|
|
829 ELSE IF( TJJ.GT.ZERO ) THEN |
|
|
830 * |
|
|
831 * 0 < abs(A(j,j)) <= SMLNUM: |
|
|
832 * |
|
|
833 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
834 * |
|
|
835 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. |
|
|
836 * |
|
|
837 REC = ( TJJ*BIGNUM ) / XJ |
|
|
838 CALL ZDSCAL( N, REC, X, 1 ) |
|
|
839 SCALE = SCALE*REC |
|
|
840 XMAX = XMAX*REC |
|
|
841 END IF |
|
|
842 X( J ) = ZLADIV( X( J ), TJJS ) |
|
|
843 ELSE |
|
|
844 * |
|
|
845 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and |
|
|
846 * scale = 0 and compute a solution to A**H *x = 0. |
|
|
847 * |
|
|
848 DO 200 I = 1, N |
|
|
849 X( I ) = ZERO |
|
|
850 200 CONTINUE |
|
|
851 X( J ) = ONE |
|
|
852 SCALE = ZERO |
|
|
853 XMAX = ZERO |
|
|
854 END IF |
|
|
855 210 CONTINUE |
|
|
856 ELSE |
|
|
857 * |
|
|
858 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot |
|
|
859 * product has already been divided by 1/A(j,j). |
|
|
860 * |
|
|
861 X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ |
|
|
862 END IF |
|
|
863 XMAX = MAX( XMAX, CABS1( X( J ) ) ) |
|
|
864 220 CONTINUE |
|
|
865 END IF |
|
|
866 SCALE = SCALE / TSCAL |
|
|
867 END IF |
|
|
868 * |
|
|
869 * Scale the column norms by 1/TSCAL for return. |
|
|
870 * |
|
|
871 IF( TSCAL.NE.ONE ) THEN |
|
|
872 CALL DSCAL( N, ONE / TSCAL, CNORM, 1 ) |
|
|
873 END IF |
|
|
874 * |
|
|
875 RETURN |
|
|
876 * |
|
|
877 * End of ZLATRS |
|
|
878 * |
|
|
879 END |
