236 SUBROUTINE slatrs( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
244 CHARACTER DIAG, NORMIN, TRANS, UPLO
249 REAL A( LDA, * ), CNORM( * ), X( * )
256 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0 )
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
267 REAL SASUM, SDOT, SLAMCH
268 EXTERNAL lsame, isamax, sasum, sdot, slamch
274 INTRINSIC abs, max, min
279 upper = lsame( uplo,
'U' )
280 notran = lsame( trans,
'N' )
281 nounit = lsame( diag,
'N' )
285 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
287 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) .AND. .NOT.
288 $ lsame( trans,
'C' ) )
THEN
290 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag,
'U' ) )
THEN
292 ELSE IF( .NOT.lsame( normin,
'Y' ) .AND. .NOT.
293 $ lsame( normin,
'N' ) )
THEN
295 ELSE IF( n.LT.0 )
THEN
297 ELSE IF( lda.LT.max( 1, n ) )
THEN
301 CALL xerbla(
'SLATRS', -info )
312 smlnum = slamch(
'Safe minimum' ) / slamch(
'Precision' )
313 bignum = one / smlnum
316 IF( lsame( normin,
'N' ) )
THEN
325 cnorm( j ) = sasum( j-1, a( 1, j ), 1 )
332 cnorm( j ) = sasum( n-j, a( j+1, j ), 1 )
341 imax = isamax( n, cnorm, 1 )
343 IF( tmax.LE.bignum )
THEN
346 tscal = one / ( smlnum*tmax )
347 CALL sscal( n, tscal, cnorm, 1 )
353 j = isamax( n, x, 1 )
370 IF( tscal.NE.one )
THEN
382 grow = one / max( xbnd, smlnum )
384 DO 30 j = jfirst, jlast, jinc
393 tjj = abs( a( j, j ) )
394 xbnd = min( xbnd, min( one, tjj )*grow )
395 IF( tjj+cnorm( j ).GE.smlnum )
THEN
399 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
414 grow = min( one, one / max( xbnd, smlnum ) )
415 DO 40 j = jfirst, jlast, jinc
424 grow = grow*( one / ( one+cnorm( j ) ) )
443 IF( tscal.NE.one )
THEN
455 grow = one / max( xbnd, smlnum )
457 DO 60 j = jfirst, jlast, jinc
466 xj = one + cnorm( j )
467 grow = min( grow, xbnd / xj )
471 tjj = abs( a( j, j ) )
473 $ xbnd = xbnd*( tjj / xj )
475 grow = min( grow, xbnd )
482 grow = min( one, one / max( xbnd, smlnum ) )
483 DO 70 j = jfirst, jlast, jinc
492 xj = one + cnorm( j )
499 IF( ( grow*tscal ).GT.smlnum )
THEN
504 CALL strsv( uplo, trans, diag, n, a, lda, x, 1 )
509 IF( xmax.GT.bignum )
THEN
514 scale = bignum / xmax
515 CALL sscal( n, scale, x, 1 )
523 DO 100 j = jfirst, jlast, jinc
529 tjjs = a( j, j )*tscal
536 IF( tjj.GT.smlnum )
THEN
540 IF( tjj.LT.one )
THEN
541 IF( xj.GT.tjj*bignum )
THEN
546 CALL sscal( n, rec, x, 1 )
551 x( j ) = x( j ) / tjjs
553 ELSE IF( tjj.GT.zero )
THEN
557 IF( xj.GT.tjj*bignum )
THEN
562 rec = ( tjj*bignum ) / xj
563 IF( cnorm( j ).GT.one )
THEN
568 rec = rec / cnorm( j )
570 CALL sscal( n, rec, x, 1 )
574 x( j ) = x( j ) / tjjs
596 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
601 CALL sscal( n, rec, x, 1 )
604 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
608 CALL sscal( n, half, x, 1 )
618 CALL saxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
620 i = isamax( j-1, x, 1 )
629 CALL saxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
631 i = j + isamax( n-j, x( j+1 ), 1 )
641 DO 140 j = jfirst, jlast, jinc
648 rec = one / max( xmax, one )
649 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
655 tjjs = a( j, j )*tscal
660 IF( tjj.GT.one )
THEN
664 rec = min( one, rec*tjj )
667 IF( rec.LT.one )
THEN
668 CALL sscal( n, rec, x, 1 )
675 IF( uscal.EQ.one )
THEN
681 sumj = sdot( j-1, a( 1, j ), 1, x, 1 )
682 ELSE IF( j.LT.n )
THEN
683 sumj = sdot( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
691 sumj = sumj + ( a( i, j )*uscal )*x( i )
693 ELSE IF( j.LT.n )
THEN
695 sumj = sumj + ( a( i, j )*uscal )*x( i )
700 IF( uscal.EQ.tscal )
THEN
705 x( j ) = x( j ) - sumj
708 tjjs = a( j, j )*tscal
718 IF( tjj.GT.smlnum )
THEN
722 IF( tjj.LT.one )
THEN
723 IF( xj.GT.tjj*bignum )
THEN
728 CALL sscal( n, rec, x, 1 )
733 x( j ) = x( j ) / tjjs
734 ELSE IF( tjj.GT.zero )
THEN
738 IF( xj.GT.tjj*bignum )
THEN
742 rec = ( tjj*bignum ) / xj
743 CALL sscal( n, rec, x, 1 )
747 x( j ) = x( j ) / tjjs
766 x( j ) = x( j ) / tjjs - sumj
768 xmax = max( xmax, abs( x( j ) ) )
771 scale = scale / tscal
776 IF( tscal.NE.one )
THEN
777 CALL sscal( n, one / tscal, cnorm, 1 )
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
SLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine strsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
STRSV