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.LE. IF( TMAXBIGNUM*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.NE. IF( TSCALONE ) 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.LE. IF( GROWSMLNUM )
404 $ GO TO 60
405*
406 TJJS = AP( IP )
407 TJJ = CABS1( TJJS )
408*
409.GE. IF( TJJSMLNUM ) 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.GE. IF( TJJ+CNORM( J )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.LE. IF( GROWSMLNUM )
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.NE. IF( TSCALONE ) 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.LE. IF( GROWSMLNUM )
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.GE. IF( TJJSMLNUM ) THEN
503*
504* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
505*
506.GT. IF( XJTJJ )
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.LE. IF( GROWSMLNUM )
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.GT. IF( ( GROW*TSCAL )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.GT. IF( XMAXBIGNUM*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.EQ. IF( TSCALONE )
578 $ GO TO 110
579 END IF
580 TJJ = CABS1( TJJS )
581.GT. IF( TJJSMLNUM ) THEN
582*
583* abs(A(j,j)) > SMLNUM:
584*
585.LT. IF( TJJONE ) THEN
586.GT. IF( XJTJJ*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.GT. ELSE IF( TJJZERO ) THEN
599*
600* 0 < abs(A(j,j)) <= SMLNUM:
601*
602.GT. IF( XJTJJ*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.GT. IF( CNORM( J )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.GT. IF( XJONE ) THEN
640 REC = ONE / XJ
641.GT. IF( CNORM( J )( 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.GT. ELSE IF( XJ*CNORM( J )( 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.GT. IF( J1 ) 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.LT. IF( JN ) 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.GT. IF( CNORM( J )( 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.GT. IF( TJJONE ) 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.LT. IF( RECONE ) 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.EQ. IF( USCALDCMPLX( 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.LT. ELSE IF( JN ) 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.LT. ELSE IF( JN ) 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.EQ. IF( USCALDCMPLX( 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.EQ. IF( TSCALONE )
764 $ GO TO 160
765 END IF
766 TJJ = CABS1( TJJS )
767.GT. IF( TJJSMLNUM ) THEN
768*
769* abs(A(j,j)) > SMLNUM:
770*
771.LT. IF( TJJONE ) THEN
772.GT. IF( XJTJJ*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.GT. ELSE IF( TJJZERO ) THEN
784*
785* 0 < abs(A(j,j)) <= SMLNUM:
786*
787.GT. IF( XJTJJ*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.GT. IF( CNORM( J )( 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.GT. IF( TJJONE ) 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.LT. IF( RECONE ) 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.EQ. IF( USCALDCMPLX( 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.LT. ELSE IF( JN ) 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.LT. ELSE IF( JN ) 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.EQ. IF( USCALDCMPLX( 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.EQ. IF( TSCALONE )
904 $ GO TO 210
905 END IF
906 TJJ = CABS1( TJJS )
907.GT. IF( TJJSMLNUM ) THEN
908*
909* abs(A(j,j)) > SMLNUM:
910*
911.LT. IF( TJJONE ) THEN
912.GT. IF( XJTJJ*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.GT. ELSE IF( TJJZERO ) THEN
924*
925* 0 < abs(A(j,j)) <= SMLNUM:
926*
927.GT. IF( XJTJJ*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.NE. IF( TSCALONE ) 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
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21