242 SUBROUTINE dlatbs( 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
253 DOUBLE PRECISION scale
256 DOUBLE PRECISION ab( ldab, * ), cnorm( * ), x( * )
262 DOUBLE PRECISION zero, half, one
263 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
266 LOGICAL notran, nounit, upper
267 INTEGER i, imax, j, jfirst, jinc, jlast, jlen, maind
268 DOUBLE PRECISION bignum, grow, rec, smlnum, sumj, tjj, tjjs,
269 $ tmax, tscal, uscal, xbnd, xj, xmax
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(
'DLATBS', -info )
321 smlnum =
dlamch(
'Safe minimum' ) /
dlamch(
'Precision' )
322 bignum = one / smlnum
325 IF(
lsame( normin,
'N' ) )
THEN
334 jlen = min( kd, j-1 )
335 cnorm( j ) =
dasum( jlen, ab( kd+1-jlen, j ), 1 )
342 jlen = min( kd, n-j )
344 cnorm( j ) =
dasum( jlen, ab( 2, j ), 1 )
355 imax =
idamax( n, cnorm, 1 )
357 IF( tmax.LE.bignum )
THEN
360 tscal = one / ( smlnum*tmax )
361 CALL
dscal( n, tscal, cnorm, 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
dtbsv( uplo, trans, diag, n, kd, ab, ldab, x, 1 )
527 IF( xmax.GT.bignum )
THEN
532 scale = bignum / xmax
533 CALL
dscal( n, scale, x, 1 )
541 DO 110 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
dscal( 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
dscal( n, rec, x, 1 )
592 x( j ) = x( j ) / tjjs
614 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
619 CALL
dscal( n, rec, x, 1 )
622 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
626 CALL
dscal( n, half, x, 1 )
637 jlen = min( kd, j-1 )
638 CALL
daxpy( jlen, -x( j )*tscal,
639 $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
643 ELSE IF( j.LT.n )
THEN
649 jlen = min( kd, n-j )
651 $ CALL
daxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
653 i = j +
idamax( n-j, x( j+1 ), 1 )
662 DO 160 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
dscal( n, rec, x, 1 )
696 IF( uscal.EQ.one )
THEN
702 jlen = min( kd, j-1 )
703 sumj =
ddot( jlen, ab( kd+1-jlen, j ), 1,
706 jlen = min( kd, n-j )
708 $ sumj =
ddot( 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
dscal( 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
dscal( 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
dscal( n, one / tscal, cnorm, 1 )