LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ slatrs()

 subroutine slatrs ( character uplo, character trans, character diag, character normin, integer n, real, dimension( lda, * ) a, integer lda, real, dimension( * ) x, real scale, real, dimension( * ) cnorm, integer info )

SLATRS solves a triangular system of equations with the scale factor set to prevent overflow.

Purpose:
``` SLATRS solves one of the triangular systems

A *x = s*b  or  A**T*x = s*b

with scaling to prevent overflow.  Here A is an upper or lower
triangular matrix, A**T denotes the transpose of A, x and b are
n-element vectors, and s is a scaling factor, usually less than
or equal to 1, chosen so that the components of x will be less than
the overflow threshold.  If the unscaled problem will not cause
overflow, the Level 2 BLAS routine STRSV is called.  If the matrix A
is singular (A(j,j) = 0 for some j), then s is set to 0 and a
non-trivial solution to A*x = 0 is returned.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangular``` [in] TRANS ``` TRANS is CHARACTER*1 Specifies the operation applied to A. = 'N': Solve A * x = s*b (No transpose) = 'T': Solve A**T* x = s*b (Transpose) = 'C': Solve A**T* x = s*b (Conjugate transpose = Transpose)``` [in] DIAG ``` DIAG is CHARACTER*1 Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangular``` [in] NORMIN ``` NORMIN is CHARACTER*1 Specifies whether CNORM has been set or not. = 'Y': CNORM contains the column norms on entry = 'N': CNORM is not set on entry. On exit, the norms will be computed and stored in CNORM.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] A ``` A is REAL array, dimension (LDA,N) The triangular matrix A. If UPLO = 'U', the leading n by n upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n by n lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. If DIAG = 'U', the diagonal elements of A are also not referenced and are assumed to be 1.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max (1,N).``` [in,out] X ``` X is REAL array, dimension (N) On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x.``` [out] SCALE ``` SCALE is REAL The scaling factor s for the triangular system A * x = s*b or A**T* x = s*b. If SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0.``` [in,out] CNORM ``` CNORM is REAL array, dimension (N) If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains the norm of the off-diagonal part of the j-th column of A. If TRANS = 'N', CNORM(j) must be greater than or equal to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be greater than or equal to the 1-norm. If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal value```
Further Details:
```  A rough bound on x is computed; if that is less than overflow, STRSV
is called, otherwise, specific code is used which checks for possible
overflow or divide-by-zero at every operation.

A columnwise scheme is used for solving A*x = b.  The basic algorithm
if A is lower triangular is

x[1:n] := b[1:n]
for j = 1, ..., n
x(j) := x(j) / A(j,j)
x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
end

Define bounds on the components of x after j iterations of the loop:
M(j) = bound on x[1:j]
G(j) = bound on x[j+1:n]
Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.

Then for iteration j+1 we have
M(j+1) <= G(j) / | A(j+1,j+1) |
G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
<= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )

where CNORM(j+1) is greater than or equal to the infinity-norm of
column j+1 of A, not counting the diagonal.  Hence

G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
1<=i<=j
and

|x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
1<=i< j

Since |x(j)| <= M(j), we use the Level 2 BLAS routine STRSV if the
reciprocal of the largest M(j), j=1,..,n, is larger than
max(underflow, 1/overflow).

The bound on x(j) is also used to determine when a step in the
columnwise method can be performed without fear of overflow.  If
the computed bound is greater than a large constant, x is scaled to
prevent overflow, but if the bound overflows, x is set to 0, x(j) to
1, and scale to 0, and a non-trivial solution to A*x = 0 is found.

Similarly, a row-wise scheme is used to solve A**T*x = b.  The basic
algorithm for A upper triangular is

for j = 1, ..., n
x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
end

We simultaneously compute two bounds
G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
M(j) = bound on x(i), 1<=i<=j

The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
Then the bound on x(j) is

M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |

<= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
1<=i<=j

and we can safely call STRSV if 1/M(n) and 1/G(n) are both greater
than max(underflow, 1/overflow).```

Definition at line 236 of file slatrs.f.

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 REAL SCALE
247* ..
248* .. Array Arguments ..
249 REAL A( LDA, * ), CNORM( * ), X( * )
250* ..
251*
252* =====================================================================
253*
254* .. Parameters ..
255 REAL ZERO, HALF, ONE
256 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0 )
257* ..
258* .. Local Scalars ..
259 LOGICAL NOTRAN, NOUNIT, UPPER
260 INTEGER I, IMAX, J, JFIRST, JINC, JLAST
261 REAL BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
262 \$ TMAX, TSCAL, USCAL, XBND, XJ, XMAX
263* ..
264* .. Local Arrays ..
265 REAL WORK(1)
266* ..
267* .. External Functions ..
268 LOGICAL LSAME
269 INTEGER ISAMAX
270 REAL SASUM, SDOT, SLAMCH, SLANGE
271 EXTERNAL lsame, isamax, sasum, sdot, slamch, slange
272* ..
273* .. External Subroutines ..
274 EXTERNAL saxpy, sscal, strsv, 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( 'SLATRS', -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 = slamch( 'Safe minimum' ) / slamch( '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 ) = sasum( 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 ) = sasum( 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 = isamax( 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.slamch('Overflow') ) THEN
354* Case 1: All entries in CNORM are valid floating-point numbers
355 tscal = one / ( smlnum*tmax )
356 CALL sscal( 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( slange( '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( slange( 'M', n-j, 1, a( j+1, j ), 1,
377 \$ work ), tmax )
378 END DO
379 END IF
380*
381 IF( tmax.LE.slamch('Overflow') ) THEN
382 tscal = one / ( smlnum*tmax )
383 DO j = 1, n
384 IF( cnorm( j ).LE.slamch('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 strsv( 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 STRSV can be used.
414*
415 j = isamax( 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 strsv( 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 sscal( n, scale, x, 1 )
578 xmax = bignum
579 END IF
580*
581 IF( notran ) THEN
582*
583* Solve A * x = b
584*
585 DO 100 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 95
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 sscal( 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 sscal( 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 95 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 sscal( 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 sscal( 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 saxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
681 \$ 1 )
682 i = isamax( 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 saxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
692 \$ x( j+1 ), 1 )
693 i = j + isamax( n-j, x( j+1 ), 1 )
694 xmax = abs( x( i ) )
695 END IF
696 END IF
697 100 CONTINUE
698*
699 ELSE
700*
701* Solve A**T * x = b
702*
703 DO 140 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 sscal( 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 SDOT to perform the dot product.
741*
742 IF( upper ) THEN
743 sumj = sdot( j-1, a( 1, j ), 1, x, 1 )
744 ELSE IF( j.LT.n ) THEN
745 sumj = sdot( 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 110 i = 1, j - 1
753 sumj = sumj + ( a( i, j )*uscal )*x( i )
754 110 CONTINUE
755 ELSE IF( j.LT.n ) THEN
756 DO 120 i = j + 1, n
757 sumj = sumj + ( a( i, j )*uscal )*x( i )
758 120 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 135
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 sscal( 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 sscal( 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 130 i = 1, n
816 x( i ) = zero
817 130 CONTINUE
818 x( j ) = one
819 scale = zero
820 xmax = zero
821 END IF
822 135 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 140 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 sscal( n, one / tscal, cnorm, 1 )
840 END IF
841*
842 RETURN
843*
844* End of SLATRS
845*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
real function sasum(n, sx, incx)
SASUM
Definition sasum.f:72
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:114
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine strsv(uplo, trans, diag, n, a, lda, x, incx)
STRSV
Definition strsv.f:149
Here is the call graph for this function:
Here is the caller graph for this function: