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