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