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