238 SUBROUTINE dlatbs( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB,
240 $ SCALE, CNORM, INFO )
247 CHARACTER DIAG, NORMIN, TRANS, UPLO
248 INTEGER INFO, KD, LDAB, N
249 DOUBLE PRECISION SCALE
252 DOUBLE PRECISION AB( LDAB, * ), CNORM( * ), X( * )
258 DOUBLE PRECISION ZERO, HALF, ONE
259 PARAMETER ( ZERO = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
262 LOGICAL NOTRAN, NOUNIT, UPPER
263 INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
264 DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
265 $ tmax, tscal, uscal, xbnd, xj, xmax
270 DOUBLE PRECISION DASUM, DDOT, DLAMCH
271 EXTERNAL lsame, idamax, dasum, ddot, dlamch
277 INTRINSIC abs, max, min
282 upper = lsame( uplo,
'U' )
283 notran = lsame( trans,
'N' )
284 nounit = lsame( diag,
'N' )
288 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
290 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) .AND. .NOT.
291 $ lsame( trans,
'C' ) )
THEN
293 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag,
'U' ) )
THEN
295 ELSE IF( .NOT.lsame( normin,
'Y' ) .AND. .NOT.
296 $ lsame( normin,
'N' ) )
THEN
298 ELSE IF( n.LT.0 )
THEN
300 ELSE IF( kd.LT.0 )
THEN
302 ELSE IF( ldab.LT.kd+1 )
THEN
306 CALL xerbla(
'DLATBS', -info )
318 smlnum = dlamch(
'Safe minimum' ) / dlamch(
'Precision' )
319 bignum = one / smlnum
321 IF( lsame( normin,
'N' ) )
THEN
330 jlen = min( kd, j-1 )
331 cnorm( j ) = dasum( jlen, ab( kd+1-jlen, j ), 1 )
338 jlen = min( kd, n-j )
340 cnorm( j ) = dasum( jlen, ab( 2, j ), 1 )
351 imax = idamax( n, cnorm, 1 )
353 IF( tmax.LE.bignum )
THEN
356 tscal = one / ( smlnum*tmax )
357 CALL dscal( n, tscal, cnorm, 1 )
363 j = idamax( n, x, 1 )
382 IF( tscal.NE.one )
THEN
394 grow = one / max( xbnd, smlnum )
396 DO 30 j = jfirst, jlast, jinc
405 tjj = abs( ab( maind, j ) )
406 xbnd = min( xbnd, min( one, tjj )*grow )
407 IF( tjj+cnorm( j ).GE.smlnum )
THEN
411 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
426 grow = min( one, one / max( xbnd, smlnum ) )
427 DO 40 j = jfirst, jlast, jinc
436 grow = grow*( one / ( one+cnorm( j ) ) )
457 IF( tscal.NE.one )
THEN
469 grow = one / max( xbnd, smlnum )
471 DO 60 j = jfirst, jlast, jinc
480 xj = one + cnorm( j )
481 grow = min( grow, xbnd / xj )
485 tjj = abs( ab( maind, j ) )
487 $ xbnd = xbnd*( tjj / xj )
489 grow = min( grow, xbnd )
496 grow = min( one, one / max( xbnd, smlnum ) )
497 DO 70 j = jfirst, jlast, jinc
506 xj = one + cnorm( j )
513 IF( ( grow*tscal ).GT.smlnum )
THEN
518 CALL dtbsv( uplo, trans, diag, n, kd, ab, ldab, x, 1 )
523 IF( xmax.GT.bignum )
THEN
528 scale = bignum / xmax
529 CALL dscal( n, scale, x, 1 )
537 DO 110 j = jfirst, jlast, jinc
543 tjjs = ab( maind, j )*tscal
550 IF( tjj.GT.smlnum )
THEN
554 IF( tjj.LT.one )
THEN
555 IF( xj.GT.tjj*bignum )
THEN
560 CALL dscal( n, rec, x, 1 )
565 x( j ) = x( j ) / tjjs
567 ELSE IF( tjj.GT.zero )
THEN
571 IF( xj.GT.tjj*bignum )
THEN
576 rec = ( tjj*bignum ) / xj
577 IF( cnorm( j ).GT.one )
THEN
582 rec = rec / cnorm( j )
584 CALL dscal( n, rec, x, 1 )
588 x( j ) = x( j ) / tjjs
610 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
615 CALL dscal( n, rec, x, 1 )
618 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
622 CALL dscal( n, half, x, 1 )
633 jlen = min( kd, j-1 )
634 CALL daxpy( jlen, -x( j )*tscal,
635 $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
636 i = idamax( j-1, x, 1 )
639 ELSE IF( j.LT.n )
THEN
645 jlen = min( kd, n-j )
647 $
CALL daxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
649 i = j + idamax( n-j, x( j+1 ), 1 )
658 DO 160 j = jfirst, jlast, jinc
665 rec = one / max( xmax, one )
666 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
672 tjjs = ab( maind, j )*tscal
677 IF( tjj.GT.one )
THEN
681 rec = min( one, rec*tjj )
684 IF( rec.LT.one )
THEN
685 CALL dscal( n, rec, x, 1 )
692 IF( uscal.EQ.one )
THEN
698 jlen = min( kd, j-1 )
699 sumj = ddot( jlen, ab( kd+1-jlen, j ), 1,
702 jlen = min( kd, n-j )
704 $ sumj = ddot( jlen, ab( 2, j ), 1, x( j+1 ),
712 jlen = min( kd, j-1 )
714 sumj = sumj + ( ab( kd+i-jlen, j )*uscal )*
718 jlen = min( kd, n-j )
720 sumj = sumj + ( ab( i+1, j )*uscal )*x( j+i )
725 IF( uscal.EQ.tscal )
THEN
730 x( j ) = x( j ) - sumj
736 tjjs = ab( maind, j )*tscal
743 IF( tjj.GT.smlnum )
THEN
747 IF( tjj.LT.one )
THEN
748 IF( xj.GT.tjj*bignum )
THEN
753 CALL dscal( n, rec, x, 1 )
758 x( j ) = x( j ) / tjjs
759 ELSE IF( tjj.GT.zero )
THEN
763 IF( xj.GT.tjj*bignum )
THEN
767 rec = ( tjj*bignum ) / xj
768 CALL dscal( n, rec, x, 1 )
772 x( j ) = x( j ) / tjjs
791 x( j ) = x( j ) / tjjs - sumj
793 xmax = max( xmax, abs( x( j ) ) )
796 scale = scale / tscal
801 IF( tscal.NE.one )
THEN
802 CALL dscal( n, one / tscal, cnorm, 1 )