242 SUBROUTINE slatbs( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
243 $ scale, cnorm, info )
251 CHARACTER DIAG, NORMIN, TRANS, UPLO
252 INTEGER INFO, KD, LDAB, N
256 REAL AB( ldab, * ), CNORM( * ), X( * )
263 parameter ( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0 )
266 LOGICAL NOTRAN, NOUNIT, UPPER
267 INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
268 REAL BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
269 $ tmax, tscal, uscal, xbnd, xj, xmax
274 REAL SASUM, SDOT, SLAMCH
275 EXTERNAL lsame, isamax, sasum, sdot, slamch
281 INTRINSIC abs, max, min
286 upper = lsame( uplo,
'U' )
287 notran = lsame( trans,
'N' )
288 nounit = lsame( diag,
'N' )
292 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
294 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) .AND. .NOT.
295 $ lsame( trans,
'C' ) )
THEN
297 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag,
'U' ) )
THEN
299 ELSE IF( .NOT.lsame( normin,
'Y' ) .AND. .NOT.
300 $ lsame( normin,
'N' ) )
THEN
302 ELSE IF( n.LT.0 )
THEN
304 ELSE IF( kd.LT.0 )
THEN
306 ELSE IF( ldab.LT.kd+1 )
THEN
310 CALL xerbla(
'SLATBS', -info )
321 smlnum = slamch(
'Safe minimum' ) / slamch(
'Precision' )
322 bignum = one / smlnum
325 IF( lsame( normin,
'N' ) )
THEN
334 jlen = min( kd, j-1 )
335 cnorm( j ) = sasum( jlen, ab( kd+1-jlen, j ), 1 )
342 jlen = min( kd, n-j )
344 cnorm( j ) = sasum( jlen, ab( 2, j ), 1 )
355 imax = isamax( n, cnorm, 1 )
357 IF( tmax.LE.bignum )
THEN
360 tscal = one / ( smlnum*tmax )
361 CALL sscal( n, tscal, cnorm, 1 )
367 j = isamax( n, x, 1 )
386 IF( tscal.NE.one )
THEN
398 grow = one / max( xbnd, smlnum )
400 DO 30 j = jfirst, jlast, jinc
409 tjj = abs( ab( maind, j ) )
410 xbnd = min( xbnd, min( one, tjj )*grow )
411 IF( tjj+cnorm( j ).GE.smlnum )
THEN
415 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
430 grow = min( one, one / max( xbnd, smlnum ) )
431 DO 40 j = jfirst, jlast, jinc
440 grow = grow*( one / ( one+cnorm( j ) ) )
461 IF( tscal.NE.one )
THEN
473 grow = one / max( xbnd, smlnum )
475 DO 60 j = jfirst, jlast, jinc
484 xj = one + cnorm( j )
485 grow = min( grow, xbnd / xj )
489 tjj = abs( ab( maind, j ) )
491 $ xbnd = xbnd*( tjj / xj )
493 grow = min( grow, xbnd )
500 grow = min( one, one / max( xbnd, smlnum ) )
501 DO 70 j = jfirst, jlast, jinc
510 xj = one + cnorm( j )
517 IF( ( grow*tscal ).GT.smlnum )
THEN
522 CALL stbsv( uplo, trans, diag, n, kd, ab, ldab, x, 1 )
527 IF( xmax.GT.bignum )
THEN
532 scale = bignum / xmax
533 CALL sscal( n, scale, x, 1 )
541 DO 100 j = jfirst, jlast, jinc
547 tjjs = ab( maind, j )*tscal
554 IF( tjj.GT.smlnum )
THEN
558 IF( tjj.LT.one )
THEN
559 IF( xj.GT.tjj*bignum )
THEN
564 CALL sscal( n, rec, x, 1 )
569 x( j ) = x( j ) / tjjs
571 ELSE IF( tjj.GT.zero )
THEN
575 IF( xj.GT.tjj*bignum )
THEN
580 rec = ( tjj*bignum ) / xj
581 IF( cnorm( j ).GT.one )
THEN
586 rec = rec / cnorm( j )
588 CALL sscal( n, rec, x, 1 )
592 x( j ) = x( j ) / tjjs
614 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
619 CALL sscal( n, rec, x, 1 )
622 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
626 CALL sscal( n, half, x, 1 )
637 jlen = min( kd, j-1 )
638 CALL saxpy( jlen, -x( j )*tscal,
639 $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
640 i = isamax( j-1, x, 1 )
643 ELSE IF( j.LT.n )
THEN
649 jlen = min( kd, n-j )
651 $
CALL saxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
653 i = j + isamax( n-j, x( j+1 ), 1 )
662 DO 140 j = jfirst, jlast, jinc
669 rec = one / max( xmax, one )
670 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
676 tjjs = ab( maind, j )*tscal
681 IF( tjj.GT.one )
THEN
685 rec = min( one, rec*tjj )
688 IF( rec.LT.one )
THEN
689 CALL sscal( n, rec, x, 1 )
696 IF( uscal.EQ.one )
THEN
702 jlen = min( kd, j-1 )
703 sumj = sdot( jlen, ab( kd+1-jlen, j ), 1,
706 jlen = min( kd, n-j )
708 $ sumj = sdot( jlen, ab( 2, j ), 1, x( j+1 ), 1 )
715 jlen = min( kd, j-1 )
717 sumj = sumj + ( ab( kd+i-jlen, j )*uscal )*
721 jlen = min( kd, n-j )
723 sumj = sumj + ( ab( i+1, j )*uscal )*x( j+i )
728 IF( uscal.EQ.tscal )
THEN
733 x( j ) = x( j ) - sumj
739 tjjs = ab( maind, j )*tscal
746 IF( tjj.GT.smlnum )
THEN
750 IF( tjj.LT.one )
THEN
751 IF( xj.GT.tjj*bignum )
THEN
756 CALL sscal( n, rec, x, 1 )
761 x( j ) = x( j ) / tjjs
762 ELSE IF( tjj.GT.zero )
THEN
766 IF( xj.GT.tjj*bignum )
THEN
770 rec = ( tjj*bignum ) / xj
771 CALL sscal( n, rec, x, 1 )
775 x( j ) = x( j ) / tjjs
794 x( j ) = x( j ) / tjjs - sumj
796 xmax = max( xmax, abs( x( j ) ) )
799 scale = scale / tscal
804 IF( tscal.NE.one )
THEN
805 CALL sscal( n, one / tscal, cnorm, 1 )
subroutine slatbs(UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
SLATBS solves a triangular banded system of equations.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine stbsv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
STBSV
subroutine sscal(N, SA, SX, INCX)
SSCAL