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