1 SUBROUTINE dlatrs( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
10 CHARACTER DIAG, NORMIN, TRANS, UPLO
12 DOUBLE PRECISION SCALE
15 DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * )
173 DOUBLE PRECISION ZERO, HALF, ONE
174 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
177 LOGICAL NOTRAN, NOUNIT, UPPER
178 INTEGER I, IMAX, J, JFIRST, JINC, JLAST
179 DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
180 $ tmax, tscal, uscal, xbnd, xj, xmax
185 DOUBLE PRECISION DASUM, DDOT, DLAMCH
186 EXTERNAL lsame, idamax, dasum, ddot, dlamch
192 INTRINSIC abs, max, min
197 upper = lsame( uplo,
'U' )
198 notran = lsame( trans,
'N' )
199 nounit = lsame( diag,
'N' )
203 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
205 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) .AND. .NOT.
206 $ lsame( trans,
'C' ) )
THEN
208 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag,
'U' ) )
THEN
210 ELSE IF( .NOT.lsame( normin,
'Y' ) .AND. .NOT.
211 $ lsame( normin,
'N' ) )
THEN
213 ELSE IF( n.LT.0 )
THEN
215 ELSE IF( lda.LT.max( 1, n ) )
THEN
219 CALL xerbla(
'DLATRS', -info )
230 smlnum = dlamch(
'Safe minimum' ) / dlamch(
'Precision' )
231 bignum = one / smlnum
234 IF( lsame( normin,
'N' ) )
THEN
243 cnorm( j ) = dasum( j-1, a( 1, j ), 1 )
250 cnorm( j ) = dasum( n-j, a( j+1, j ), 1 )
259 imax = idamax( n, cnorm, 1 )
261 IF( tmax.LE.bignum )
THEN
264 tscal = one / ( smlnum*tmax )
265 CALL dscal( n, tscal, cnorm, 1 )
271 j = idamax( n, x, 1 )
288 IF( tscal.NE.one )
THEN
300 grow = one / max( xbnd, smlnum )
302 DO 30 j = jfirst, jlast, jinc
311 tjj = abs( a( j, j ) )
312 xbnd = min( xbnd, min( one, tjj )*grow )
313 IF( tjj+cnorm( j ).GE.smlnum )
THEN
317 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
332 grow = min( one, one / max( xbnd, smlnum ) )
333 DO 40 j = jfirst, jlast, jinc
342 grow = grow*( one / ( one+cnorm( j ) ) )
361 IF( tscal.NE.one )
THEN
373 grow = one / max( xbnd, smlnum )
375 DO 60 j = jfirst, jlast, jinc
384 xj = one + cnorm( j )
385 grow = min( grow, xbnd / xj )
389 tjj = abs( a( j, j ) )
391 $ xbnd = xbnd*( tjj / xj )
393 grow = min( grow, xbnd )
400 grow = min( one, one / max( xbnd, smlnum ) )
401 DO 70 j = jfirst, jlast, jinc
410 xj = one + cnorm( j )
417 IF( ( grow*tscal ).GT.smlnum )
THEN
422 CALL dtrsv( uplo, trans, diag, n, a, lda, x, 1 )
427 IF( xmax.GT.bignum )
THEN
432 scale = bignum / xmax
433 CALL dscal( n, scale, x, 1 )
441 DO 110 j = jfirst, jlast, jinc
447 tjjs = a( j, j )*tscal
454 IF( tjj.GT.smlnum )
THEN
458 IF( tjj.LT.one )
THEN
459 IF( xj.GT.tjj*bignum )
THEN
464 CALL dscal( n, rec, x, 1 )
469 x( j ) = x( j ) / tjjs
471 ELSE IF( tjj.GT.zero )
THEN
475 IF( xj.GT.tjj*bignum )
THEN
480 rec = ( tjj*bignum ) / xj
481 IF( cnorm( j ).GT.one )
THEN
486 rec = rec / cnorm( j )
488 CALL dscal( n, rec, x, 1 )
492 x( j ) = x( j ) / tjjs
514 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
519 CALL dscal( n, rec, x, 1 )
522 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
526 CALL dscal( n, half, x, 1 )
536 CALL daxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
538 i = idamax( j-1, x, 1 )
547 CALL daxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
549 i = j + idamax( n-j, x( j+1 ), 1 )
559 DO 160 j = jfirst, jlast, jinc
566 rec = one / max( xmax, one )
567 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
573 tjjs = a( j, j )*tscal
578 IF( tjj.GT.one )
THEN
582 rec = min( one, rec*tjj )
585 IF( rec.LT.one )
THEN
586 CALL dscal( n, rec, x, 1 )
593 IF( uscal.EQ.one )
THEN
599 sumj = ddot( j-1, a( 1, j ), 1, x, 1 )
600 ELSE IF( j.LT.n )
THEN
601 sumj = ddot( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
609 sumj = sumj + ( a( i, j )*uscal )*x( i )
611 ELSE IF( j.LT.n )
THEN
613 sumj = sumj + ( a( i, j )*uscal )*x( i )
618 IF( uscal.EQ.tscal )
THEN
623 x( j ) = x( j ) - sumj
626 tjjs = a( j, j )*tscal
636 IF( tjj.GT.smlnum )
THEN
640 IF( tjj.LT.one )
THEN
641 IF( xj.GT.tjj*bignum )
THEN
646 CALL dscal( n, rec, x, 1 )
651 x( j ) = x( j ) / tjjs
652 ELSE IF( tjj.GT.zero )
THEN
656 IF( xj.GT.tjj*bignum )
THEN
660 rec = ( tjj*bignum ) / xj
661 CALL dscal( n, rec, x, 1 )
665 x( j ) = x( j ) / tjjs
684 x( j ) = x( j ) / tjjs - sumj
686 xmax = max( xmax, abs( x( j ) ) )
689 scale = scale / tscal
694 IF( tscal.NE.one )
THEN
695 CALL dscal( n, one / tscal, cnorm, 1 )
subroutine daxpy(n, da, dx, incx, dy, incy)
subroutine dlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
subroutine dscal(n, da, dx, incx)
subroutine dtrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
subroutine xerbla(SRNAME, INFO)