238 SUBROUTINE dlatrs( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
247 CHARACTER DIAG, NORMIN, TRANS, UPLO
249 DOUBLE PRECISION SCALE
252 DOUBLE PRECISION A( lda, * ), 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
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( lda.LT.max( 1, n ) )
THEN
304 CALL xerbla(
'DLATRS', -info )
315 smlnum = dlamch(
'Safe minimum' ) / dlamch(
'Precision' )
316 bignum = one / smlnum
319 IF( lsame( normin,
'N' ) )
THEN
328 cnorm( j ) = dasum( j-1, a( 1, j ), 1 )
335 cnorm( j ) = dasum( n-j, a( j+1, j ), 1 )
344 imax = idamax( n, cnorm, 1 )
346 IF( tmax.LE.bignum )
THEN
349 tscal = one / ( smlnum*tmax )
350 CALL dscal( n, tscal, cnorm, 1 )
356 j = idamax( n, x, 1 )
373 IF( tscal.NE.one )
THEN
385 grow = one / max( xbnd, smlnum )
387 DO 30 j = jfirst, jlast, jinc
396 tjj = abs( a( j, j ) )
397 xbnd = min( xbnd, min( one, tjj )*grow )
398 IF( tjj+cnorm( j ).GE.smlnum )
THEN
402 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
417 grow = min( one, one / max( xbnd, smlnum ) )
418 DO 40 j = jfirst, jlast, jinc
427 grow = grow*( one / ( one+cnorm( j ) ) )
446 IF( tscal.NE.one )
THEN
458 grow = one / max( xbnd, smlnum )
460 DO 60 j = jfirst, jlast, jinc
469 xj = one + cnorm( j )
470 grow = min( grow, xbnd / xj )
474 tjj = abs( a( j, j ) )
476 $ xbnd = xbnd*( tjj / xj )
478 grow = min( grow, xbnd )
485 grow = min( one, one / max( xbnd, smlnum ) )
486 DO 70 j = jfirst, jlast, jinc
495 xj = one + cnorm( j )
502 IF( ( grow*tscal ).GT.smlnum )
THEN
507 CALL dtrsv( uplo, trans, diag, n, a, lda, x, 1 )
512 IF( xmax.GT.bignum )
THEN
517 scale = bignum / xmax
518 CALL dscal( n, scale, x, 1 )
526 DO 110 j = jfirst, jlast, jinc
532 tjjs = a( j, j )*tscal
539 IF( tjj.GT.smlnum )
THEN
543 IF( tjj.LT.one )
THEN
544 IF( xj.GT.tjj*bignum )
THEN
549 CALL dscal( n, rec, x, 1 )
554 x( j ) = x( j ) / tjjs
556 ELSE IF( tjj.GT.zero )
THEN
560 IF( xj.GT.tjj*bignum )
THEN
565 rec = ( tjj*bignum ) / xj
566 IF( cnorm( j ).GT.one )
THEN
571 rec = rec / cnorm( j )
573 CALL dscal( n, rec, x, 1 )
577 x( j ) = x( j ) / tjjs
599 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
604 CALL dscal( n, rec, x, 1 )
607 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
611 CALL dscal( n, half, x, 1 )
621 CALL daxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
623 i = idamax( j-1, x, 1 )
632 CALL daxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
634 i = j + idamax( n-j, x( j+1 ), 1 )
644 DO 160 j = jfirst, jlast, jinc
651 rec = one / max( xmax, one )
652 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
658 tjjs = a( j, j )*tscal
663 IF( tjj.GT.one )
THEN
667 rec = min( one, rec*tjj )
670 IF( rec.LT.one )
THEN
671 CALL dscal( n, rec, x, 1 )
678 IF( uscal.EQ.one )
THEN
684 sumj = ddot( j-1, a( 1, j ), 1, x, 1 )
685 ELSE IF( j.LT.n )
THEN
686 sumj = ddot( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
694 sumj = sumj + ( a( i, j )*uscal )*x( i )
696 ELSE IF( j.LT.n )
THEN
698 sumj = sumj + ( a( i, j )*uscal )*x( i )
703 IF( uscal.EQ.tscal )
THEN
708 x( j ) = x( j ) - sumj
711 tjjs = a( j, j )*tscal
721 IF( tjj.GT.smlnum )
THEN
725 IF( tjj.LT.one )
THEN
726 IF( xj.GT.tjj*bignum )
THEN
731 CALL dscal( n, rec, x, 1 )
736 x( j ) = x( j ) / tjjs
737 ELSE IF( tjj.GT.zero )
THEN
741 IF( xj.GT.tjj*bignum )
THEN
745 rec = ( tjj*bignum ) / xj
746 CALL dscal( n, rec, x, 1 )
750 x( j ) = x( j ) / tjjs
769 x( j ) = x( j ) / tjjs - sumj
771 xmax = max( xmax, abs( x( j ) ) )
774 scale = scale / tscal
779 IF( tscal.NE.one )
THEN
780 CALL dscal( n, one / tscal, cnorm, 1 )
subroutine dlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
DLATRS solves a triangular system of equations with the scale factor set to prevent overflow...
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dtrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRSV