236 SUBROUTINE dlatrs( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
244 CHARACTER DIAG, NORMIN, TRANS, UPLO
246 DOUBLE PRECISION SCALE
249 DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * )
255 DOUBLE PRECISION ZERO, HALF, ONE
256 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
259 LOGICAL NOTRAN, NOUNIT, UPPER
260 INTEGER I, IMAX, J, JFIRST, JINC, JLAST
261 DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
262 $ tmax, tscal, uscal, xbnd, xj, xmax
267 DOUBLE PRECISION DASUM, DDOT, DLAMCH
268 EXTERNAL lsame, idamax, dasum, ddot, dlamch
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(
'DLATRS', -info )
312 smlnum = dlamch(
'Safe minimum' ) / dlamch(
'Precision' )
313 bignum = one / smlnum
316 IF( lsame( normin,
'N' ) )
THEN
325 cnorm( j ) = dasum( j-1, a( 1, j ), 1 )
332 cnorm( j ) = dasum( n-j, a( j+1, j ), 1 )
341 imax = idamax( n, cnorm, 1 )
343 IF( tmax.LE.bignum )
THEN
346 tscal = one / ( smlnum*tmax )
347 CALL dscal( n, tscal, cnorm, 1 )
353 j = idamax( 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 dtrsv( uplo, trans, diag, n, a, lda, x, 1 )
509 IF( xmax.GT.bignum )
THEN
514 scale = bignum / xmax
515 CALL dscal( n, scale, x, 1 )
523 DO 110 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 dscal( 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 dscal( n, rec, x, 1 )
574 x( j ) = x( j ) / tjjs
596 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
601 CALL dscal( n, rec, x, 1 )
604 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
608 CALL dscal( n, half, x, 1 )
618 CALL daxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
620 i = idamax( j-1, x, 1 )
629 CALL daxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
631 i = j + idamax( n-j, x( j+1 ), 1 )
641 DO 160 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 dscal( n, rec, x, 1 )
675 IF( uscal.EQ.one )
THEN
681 sumj = ddot( j-1, a( 1, j ), 1, x, 1 )
682 ELSE IF( j.LT.n )
THEN
683 sumj = ddot( 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 dscal( 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 dscal( 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 dscal( n, one / tscal, cnorm, 1 )
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dtrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRSV
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.