LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clatps.f
Go to the documentation of this file.
1*> \brief \b CLATPS 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 CLATPS + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clatps.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clatps.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clatps.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLATPS( 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* REAL SCALE
28* ..
29* .. Array Arguments ..
30* REAL CNORM( * )
31* COMPLEX AP( * ), X( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> CLATPS 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 CTPSV 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 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 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 REAL
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 REAL 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 latps
155*
156*> \par Further Details:
157* =====================
158*>
159*> \verbatim
160*>
161*> A rough bound on x is computed; if that is less than overflow, CTPSV
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 CTPSV 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 CTPSV if 1/M(n) and 1/G(n) are both greater
225*> than max(underflow, 1/overflow).
226*> \endverbatim
227*>
228* =====================================================================
229 SUBROUTINE clatps( 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 REAL SCALE
240* ..
241* .. Array Arguments ..
242 REAL CNORM( * )
243 COMPLEX AP( * ), X( * )
244* ..
245*
246* =====================================================================
247*
248* .. Parameters ..
249 REAL ZERO, HALF, ONE, TWO
250 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
251 $ two = 2.0e+0 )
252* ..
253* .. Local Scalars ..
254 LOGICAL NOTRAN, NOUNIT, UPPER
255 INTEGER I, IMAX, IP, J, JFIRST, JINC, JLAST, JLEN
256 REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
257 $ xbnd, xj, xmax
258 COMPLEX CSUMJ, TJJS, USCAL, ZDUM
259* ..
260* .. External Functions ..
261 LOGICAL LSAME
262 INTEGER ICAMAX, ISAMAX
263 REAL SCASUM, SLAMCH
264 COMPLEX CDOTC, CDOTU, CLADIV
265 EXTERNAL lsame, icamax, isamax, scasum, slamch, cdotc,
266 $ cdotu, cladiv
267* ..
268* .. External Subroutines ..
269 EXTERNAL caxpy, csscal, ctpsv, sscal, xerbla
270* ..
271* .. Intrinsic Functions ..
272 INTRINSIC abs, aimag, cmplx, conjg, max, min, real
273* ..
274* .. Statement Functions ..
275 REAL CABS1, CABS2
276* ..
277* .. Statement Function definitions ..
278 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
279 cabs2( zdum ) = abs( real( zdum ) / 2. ) +
280 $ abs( aimag( zdum ) / 2. )
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( 'CLATPS', -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 = slamch( 'Safe minimum' )
317 bignum = one / smlnum
318 smlnum = smlnum / slamch( 'Precision' )
319 bignum = one / smlnum
320 scale = one
321*
322 IF( lsame( normin, 'N' ) ) THEN
323*
324* Compute the 1-norm of each column, not including the diagonal.
325*
326 IF( upper ) THEN
327*
328* A is upper triangular.
329*
330 ip = 1
331 DO 10 j = 1, n
332 cnorm( j ) = scasum( j-1, ap( ip ), 1 )
333 ip = ip + j
334 10 CONTINUE
335 ELSE
336*
337* A is lower triangular.
338*
339 ip = 1
340 DO 20 j = 1, n - 1
341 cnorm( j ) = scasum( n-j, ap( ip+1 ), 1 )
342 ip = ip + n - j + 1
343 20 CONTINUE
344 cnorm( n ) = zero
345 END IF
346 END IF
347*
348* Scale the column norms by TSCAL if the maximum element in CNORM is
349* greater than BIGNUM/2.
350*
351 imax = isamax( n, cnorm, 1 )
352 tmax = cnorm( imax )
353 IF( tmax.LE.bignum*half ) THEN
354 tscal = one
355 ELSE
356 tscal = half / ( smlnum*tmax )
357 CALL sscal( n, tscal, cnorm, 1 )
358 END IF
359*
360* Compute a bound on the computed solution vector to see if the
361* Level 2 BLAS routine CTPSV can be used.
362*
363 xmax = zero
364 DO 30 j = 1, n
365 xmax = max( xmax, cabs2( x( j ) ) )
366 30 CONTINUE
367 xbnd = xmax
368 IF( notran ) THEN
369*
370* Compute the growth in A * x = b.
371*
372 IF( upper ) THEN
373 jfirst = n
374 jlast = 1
375 jinc = -1
376 ELSE
377 jfirst = 1
378 jlast = n
379 jinc = 1
380 END IF
381*
382 IF( tscal.NE.one ) THEN
383 grow = zero
384 GO TO 60
385 END IF
386*
387 IF( nounit ) THEN
388*
389* A is non-unit triangular.
390*
391* Compute GROW = 1/G(j) and XBND = 1/M(j).
392* Initially, G(0) = max{x(i), i=1,...,n}.
393*
394 grow = half / max( xbnd, smlnum )
395 xbnd = grow
396 ip = jfirst*( jfirst+1 ) / 2
397 jlen = n
398 DO 40 j = jfirst, jlast, jinc
399*
400* Exit the loop if the growth factor is too small.
401*
402 IF( grow.LE.smlnum )
403 $ GO TO 60
404*
405 tjjs = ap( ip )
406 tjj = cabs1( tjjs )
407*
408 IF( tjj.GE.smlnum ) THEN
409*
410* M(j) = G(j-1) / abs(A(j,j))
411*
412 xbnd = min( xbnd, min( one, tjj )*grow )
413 ELSE
414*
415* M(j) could overflow, set XBND to 0.
416*
417 xbnd = zero
418 END IF
419*
420 IF( tjj+cnorm( j ).GE.smlnum ) THEN
421*
422* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
423*
424 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
425 ELSE
426*
427* G(j) could overflow, set GROW to 0.
428*
429 grow = zero
430 END IF
431 ip = ip + jinc*jlen
432 jlen = jlen - 1
433 40 CONTINUE
434 grow = xbnd
435 ELSE
436*
437* A is unit triangular.
438*
439* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
440*
441 grow = min( one, half / max( xbnd, smlnum ) )
442 DO 50 j = jfirst, jlast, jinc
443*
444* Exit the loop if the growth factor is too small.
445*
446 IF( grow.LE.smlnum )
447 $ GO TO 60
448*
449* G(j) = G(j-1)*( 1 + CNORM(j) )
450*
451 grow = grow*( one / ( one+cnorm( j ) ) )
452 50 CONTINUE
453 END IF
454 60 CONTINUE
455*
456 ELSE
457*
458* Compute the growth in A**T * x = b or A**H * x = b.
459*
460 IF( upper ) THEN
461 jfirst = 1
462 jlast = n
463 jinc = 1
464 ELSE
465 jfirst = n
466 jlast = 1
467 jinc = -1
468 END IF
469*
470 IF( tscal.NE.one ) THEN
471 grow = zero
472 GO TO 90
473 END IF
474*
475 IF( nounit ) THEN
476*
477* A is non-unit triangular.
478*
479* Compute GROW = 1/G(j) and XBND = 1/M(j).
480* Initially, M(0) = max{x(i), i=1,...,n}.
481*
482 grow = half / max( xbnd, smlnum )
483 xbnd = grow
484 ip = jfirst*( jfirst+1 ) / 2
485 jlen = 1
486 DO 70 j = jfirst, jlast, jinc
487*
488* Exit the loop if the growth factor is too small.
489*
490 IF( grow.LE.smlnum )
491 $ GO TO 90
492*
493* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
494*
495 xj = one + cnorm( j )
496 grow = min( grow, xbnd / xj )
497*
498 tjjs = ap( ip )
499 tjj = cabs1( tjjs )
500*
501 IF( tjj.GE.smlnum ) THEN
502*
503* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
504*
505 IF( xj.GT.tjj )
506 $ xbnd = xbnd*( tjj / xj )
507 ELSE
508*
509* M(j) could overflow, set XBND to 0.
510*
511 xbnd = zero
512 END IF
513 jlen = jlen + 1
514 ip = ip + jinc*jlen
515 70 CONTINUE
516 grow = min( grow, xbnd )
517 ELSE
518*
519* A is unit triangular.
520*
521* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
522*
523 grow = min( one, half / max( xbnd, smlnum ) )
524 DO 80 j = jfirst, jlast, jinc
525*
526* Exit the loop if the growth factor is too small.
527*
528 IF( grow.LE.smlnum )
529 $ GO TO 90
530*
531* G(j) = ( 1 + CNORM(j) )*G(j-1)
532*
533 xj = one + cnorm( j )
534 grow = grow / xj
535 80 CONTINUE
536 END IF
537 90 CONTINUE
538 END IF
539*
540 IF( ( grow*tscal ).GT.smlnum ) THEN
541*
542* Use the Level 2 BLAS solve if the reciprocal of the bound on
543* elements of X is not too small.
544*
545 CALL ctpsv( uplo, trans, diag, n, ap, x, 1 )
546 ELSE
547*
548* Use a Level 1 BLAS solve, scaling intermediate results.
549*
550 IF( xmax.GT.bignum*half ) THEN
551*
552* Scale X so that its components are less than or equal to
553* BIGNUM in absolute value.
554*
555 scale = ( bignum*half ) / xmax
556 CALL csscal( n, scale, x, 1 )
557 xmax = bignum
558 ELSE
559 xmax = xmax*two
560 END IF
561*
562 IF( notran ) THEN
563*
564* Solve A * x = b
565*
566 ip = jfirst*( jfirst+1 ) / 2
567 DO 110 j = jfirst, jlast, jinc
568*
569* Compute x(j) = b(j) / A(j,j), scaling x if necessary.
570*
571 xj = cabs1( x( j ) )
572 IF( nounit ) THEN
573 tjjs = ap( ip )*tscal
574 ELSE
575 tjjs = tscal
576 IF( tscal.EQ.one )
577 $ GO TO 105
578 END IF
579 tjj = cabs1( tjjs )
580 IF( tjj.GT.smlnum ) THEN
581*
582* abs(A(j,j)) > SMLNUM:
583*
584 IF( tjj.LT.one ) THEN
585 IF( xj.GT.tjj*bignum ) THEN
586*
587* Scale x by 1/b(j).
588*
589 rec = one / xj
590 CALL csscal( n, rec, x, 1 )
591 scale = scale*rec
592 xmax = xmax*rec
593 END IF
594 END IF
595 x( j ) = cladiv( x( j ), tjjs )
596 xj = cabs1( x( j ) )
597 ELSE IF( tjj.GT.zero ) THEN
598*
599* 0 < abs(A(j,j)) <= SMLNUM:
600*
601 IF( xj.GT.tjj*bignum ) THEN
602*
603* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
604* to avoid overflow when dividing by A(j,j).
605*
606 rec = ( tjj*bignum ) / xj
607 IF( cnorm( j ).GT.one ) THEN
608*
609* Scale by 1/CNORM(j) to avoid overflow when
610* multiplying x(j) times column j.
611*
612 rec = rec / cnorm( j )
613 END IF
614 CALL csscal( n, rec, x, 1 )
615 scale = scale*rec
616 xmax = xmax*rec
617 END IF
618 x( j ) = cladiv( x( j ), tjjs )
619 xj = cabs1( x( j ) )
620 ELSE
621*
622* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
623* scale = 0, and compute a solution to A*x = 0.
624*
625 DO 100 i = 1, n
626 x( i ) = zero
627 100 CONTINUE
628 x( j ) = one
629 xj = one
630 scale = zero
631 xmax = zero
632 END IF
633 105 CONTINUE
634*
635* Scale x if necessary to avoid overflow when adding a
636* multiple of column j of A.
637*
638 IF( xj.GT.one ) THEN
639 rec = one / xj
640 IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
641*
642* Scale x by 1/(2*abs(x(j))).
643*
644 rec = rec*half
645 CALL csscal( n, rec, x, 1 )
646 scale = scale*rec
647 END IF
648 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
649*
650* Scale x by 1/2.
651*
652 CALL csscal( n, half, x, 1 )
653 scale = scale*half
654 END IF
655*
656 IF( upper ) THEN
657 IF( j.GT.1 ) THEN
658*
659* Compute the update
660* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
661*
662 CALL caxpy( j-1, -x( j )*tscal, ap( ip-j+1 ), 1, x,
663 $ 1 )
664 i = icamax( 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 caxpy( n-j, -x( j )*tscal, ap( ip+1 ), 1,
675 $ x( j+1 ), 1 )
676 i = j + icamax( n-j, x( j+1 ), 1 )
677 xmax = cabs1( x( i ) )
678 END IF
679 ip = ip + n - j + 1
680 END IF
681 110 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 150 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 = cladiv( uscal, tjjs )
714 END IF
715 IF( rec.LT.one ) THEN
716 CALL csscal( 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.cmplx( one ) ) THEN
724*
725* If the scaling needed for A in the dot product is 1,
726* call CDOTU to perform the dot product.
727*
728 IF( upper ) THEN
729 csumj = cdotu( j-1, ap( ip-j+1 ), 1, x, 1 )
730 ELSE IF( j.LT.n ) THEN
731 csumj = cdotu( 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 120 i = 1, j - 1
739 csumj = csumj + ( ap( ip-j+i )*uscal )*x( i )
740 120 CONTINUE
741 ELSE IF( j.LT.n ) THEN
742 DO 130 i = 1, n - j
743 csumj = csumj + ( ap( ip+i )*uscal )*x( j+i )
744 130 CONTINUE
745 END IF
746 END IF
747*
748 IF( uscal.EQ.cmplx( 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 145
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 csscal( n, rec, x, 1 )
777 scale = scale*rec
778 xmax = xmax*rec
779 END IF
780 END IF
781 x( j ) = cladiv( 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 csscal( n, rec, x, 1 )
792 scale = scale*rec
793 xmax = xmax*rec
794 END IF
795 x( j ) = cladiv( 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 140 i = 1, n
802 x( i ) = zero
803 140 CONTINUE
804 x( j ) = one
805 scale = zero
806 xmax = zero
807 END IF
808 145 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 ) = cladiv( x( j ), tjjs ) - csumj
815 END IF
816 xmax = max( xmax, cabs1( x( j ) ) )
817 jlen = jlen + 1
818 ip = ip + jinc*jlen
819 150 CONTINUE
820*
821 ELSE
822*
823* Solve A**H * x = b
824*
825 ip = jfirst*( jfirst+1 ) / 2
826 jlen = 1
827 DO 190 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 = conjg( 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 = cladiv( uscal, tjjs )
852 END IF
853 IF( rec.LT.one ) THEN
854 CALL csscal( 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.cmplx( one ) ) THEN
862*
863* If the scaling needed for A in the dot product is 1,
864* call CDOTC to perform the dot product.
865*
866 IF( upper ) THEN
867 csumj = cdotc( j-1, ap( ip-j+1 ), 1, x, 1 )
868 ELSE IF( j.LT.n ) THEN
869 csumj = cdotc( 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 160 i = 1, j - 1
877 csumj = csumj + ( conjg( ap( ip-j+i ) )*uscal )*
878 $ x( i )
879 160 CONTINUE
880 ELSE IF( j.LT.n ) THEN
881 DO 170 i = 1, n - j
882 csumj = csumj + ( conjg( ap( ip+i ) )*uscal )*
883 $ x( j+i )
884 170 CONTINUE
885 END IF
886 END IF
887*
888 IF( uscal.EQ.cmplx( 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 = conjg( ap( ip ) )*tscal
900 ELSE
901 tjjs = tscal
902 IF( tscal.EQ.one )
903 $ GO TO 185
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 csscal( n, rec, x, 1 )
917 scale = scale*rec
918 xmax = xmax*rec
919 END IF
920 END IF
921 x( j ) = cladiv( 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 csscal( n, rec, x, 1 )
932 scale = scale*rec
933 xmax = xmax*rec
934 END IF
935 x( j ) = cladiv( 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 180 i = 1, n
942 x( i ) = zero
943 180 CONTINUE
944 x( j ) = one
945 scale = zero
946 xmax = zero
947 END IF
948 185 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 ) = cladiv( x( j ), tjjs ) - csumj
955 END IF
956 xmax = max( xmax, cabs1( x( j ) ) )
957 jlen = jlen + 1
958 ip = ip + jinc*jlen
959 190 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 sscal( n, one / tscal, cnorm, 1 )
968 END IF
969*
970 RETURN
971*
972* End of CLATPS
973*
974 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
clatps
subroutine clatps(uplo, trans, diag, normin, n, ap, x, scale, cnorm, info)
CLATPS solves a triangular system of equations with the matrix held in packed storage.
Definition clatps.f:231
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
ctpsv
subroutine ctpsv(uplo, trans, diag, n, ap, x, incx)
CTPSV
Definition ctpsv.f:144