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