237 SUBROUTINE zlatrs( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
245 CHARACTER DIAG, NORMIN, TRANS, UPLO
247 DOUBLE PRECISION SCALE
250 DOUBLE PRECISION CNORM( * )
251 COMPLEX*16 A( LDA, * ), X( * )
257 DOUBLE PRECISION ZERO, HALF, ONE, TWO
258 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, TJJ, TMAX, TSCAL,
266 COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
270 INTEGER IDAMAX, IZAMAX
271 DOUBLE PRECISION DLAMCH, DZASUM
272 COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
273 EXTERNAL lsame, idamax, izamax, dlamch, dzasum, zdotc,
280 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
283 DOUBLE PRECISION CABS1, CABS2
286 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
287 cabs2( zdum ) = abs( dble( zdum ) / 2.d0 ) +
288 $ abs( dimag( zdum ) / 2.d0 )
293 upper = lsame( uplo,
'U' )
294 notran = lsame( trans,
'N' )
295 nounit = lsame( diag,
'N' )
299 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
301 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) .AND. .NOT.
302 $ lsame( trans,
'C' ) )
THEN
304 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag,
'U' ) )
THEN
306 ELSE IF( .NOT.lsame( normin,
'Y' ) .AND. .NOT.
307 $ lsame( normin,
'N' ) )
THEN
309 ELSE IF( n.LT.0 )
THEN
311 ELSE IF( lda.LT.max( 1, n ) )
THEN
315 CALL xerbla(
'ZLATRS', -info )
326 smlnum = dlamch(
'Safe minimum' )
327 bignum = one / smlnum
328 CALL dlabad( smlnum, bignum )
329 smlnum = smlnum / dlamch(
'Precision' )
330 bignum = one / smlnum
333 IF( lsame( normin,
'N' ) )
THEN
342 cnorm( j ) = dzasum( j-1, a( 1, j ), 1 )
349 cnorm( j ) = dzasum( n-j, a( j+1, j ), 1 )
358 imax = idamax( n, cnorm, 1 )
360 IF( tmax.LE.bignum*half )
THEN
363 tscal = half / ( smlnum*tmax )
364 CALL dscal( n, tscal, cnorm, 1 )
372 xmax = max( xmax, cabs2( x( j ) ) )
390 IF( tscal.NE.one )
THEN
402 grow = half / max( xbnd, smlnum )
404 DO 40 j = jfirst, jlast, jinc
414 IF( tjj.GE.smlnum )
THEN
418 xbnd = min( xbnd, min( one, tjj )*grow )
426 IF( tjj+cnorm( j ).GE.smlnum )
THEN
430 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
445 grow = min( one, half / max( xbnd, smlnum ) )
446 DO 50 j = jfirst, jlast, jinc
455 grow = grow*( one / ( one+cnorm( j ) ) )
474 IF( tscal.NE.one )
THEN
486 grow = half / max( xbnd, smlnum )
488 DO 70 j = jfirst, jlast, jinc
497 xj = one + cnorm( j )
498 grow = min( grow, xbnd / xj )
503 IF( tjj.GE.smlnum )
THEN
508 $ xbnd = xbnd*( tjj / xj )
516 grow = min( grow, xbnd )
523 grow = min( one, half / max( xbnd, smlnum ) )
524 DO 80 j = jfirst, jlast, jinc
533 xj = one + cnorm( j )
540 IF( ( grow*tscal ).GT.smlnum )
THEN
545 CALL ztrsv( uplo, trans, diag, n, a, lda, x, 1 )
550 IF( xmax.GT.bignum*half )
THEN
555 scale = ( bignum*half ) / xmax
556 CALL zdscal( n, scale, x, 1 )
566 DO 120 j = jfirst, jlast, jinc
572 tjjs = a( j, j )*tscal
579 IF( tjj.GT.smlnum )
THEN
583 IF( tjj.LT.one )
THEN
584 IF( xj.GT.tjj*bignum )
THEN
589 CALL zdscal( n, rec, x, 1 )
594 x( j ) = zladiv( x( j ), tjjs )
596 ELSE IF( tjj.GT.zero )
THEN
600 IF( xj.GT.tjj*bignum )
THEN
605 rec = ( tjj*bignum ) / xj
606 IF( cnorm( j ).GT.one )
THEN
611 rec = rec / cnorm( j )
613 CALL zdscal( n, rec, x, 1 )
617 x( j ) = zladiv( x( j ), tjjs )
639 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
644 CALL zdscal( n, rec, x, 1 )
647 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
651 CALL zdscal( n, half, x, 1 )
661 CALL zaxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
663 i = izamax( j-1, x, 1 )
664 xmax = cabs1( x( i ) )
672 CALL zaxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
674 i = j + izamax( n-j, x( j+1 ), 1 )
675 xmax = cabs1( x( i ) )
680 ELSE IF( lsame( trans,
'T' ) )
THEN
684 DO 170 j = jfirst, jlast, jinc
691 rec = one / max( xmax, one )
692 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
698 tjjs = a( j, j )*tscal
703 IF( tjj.GT.one )
THEN
707 rec = min( one, rec*tjj )
708 uscal = zladiv( uscal, tjjs )
710 IF( rec.LT.one )
THEN
711 CALL zdscal( n, rec, x, 1 )
718 IF( uscal.EQ.dcmplx( one ) )
THEN
724 csumj = zdotu( j-1, a( 1, j ), 1, x, 1 )
725 ELSE IF( j.LT.n )
THEN
726 csumj = zdotu( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
734 csumj = csumj + ( a( i, j )*uscal )*x( i )
736 ELSE IF( j.LT.n )
THEN
738 csumj = csumj + ( a( i, j )*uscal )*x( i )
743 IF( uscal.EQ.dcmplx( tscal ) )
THEN
748 x( j ) = x( j ) - csumj
751 tjjs = a( j, j )*tscal
761 IF( tjj.GT.smlnum )
THEN
765 IF( tjj.LT.one )
THEN
766 IF( xj.GT.tjj*bignum )
THEN
771 CALL zdscal( n, rec, x, 1 )
776 x( j ) = zladiv( x( j ), tjjs )
777 ELSE IF( tjj.GT.zero )
THEN
781 IF( xj.GT.tjj*bignum )
THEN
785 rec = ( tjj*bignum ) / xj
786 CALL zdscal( n, rec, x, 1 )
790 x( j ) = zladiv( x( j ), tjjs )
809 x( j ) = zladiv( x( j ), tjjs ) - csumj
811 xmax = max( xmax, cabs1( x( j ) ) )
818 DO 220 j = jfirst, jlast, jinc
825 rec = one / max( xmax, one )
826 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
832 tjjs = dconjg( a( j, j ) )*tscal
837 IF( tjj.GT.one )
THEN
841 rec = min( one, rec*tjj )
842 uscal = zladiv( uscal, tjjs )
844 IF( rec.LT.one )
THEN
845 CALL zdscal( n, rec, x, 1 )
852 IF( uscal.EQ.dcmplx( one ) )
THEN
858 csumj = zdotc( j-1, a( 1, j ), 1, x, 1 )
859 ELSE IF( j.LT.n )
THEN
860 csumj = zdotc( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
868 csumj = csumj + ( dconjg( a( i, j ) )*uscal )*
871 ELSE IF( j.LT.n )
THEN
873 csumj = csumj + ( dconjg( a( i, j ) )*uscal )*
879 IF( uscal.EQ.dcmplx( tscal ) )
THEN
884 x( j ) = x( j ) - csumj
887 tjjs = dconjg( a( j, j ) )*tscal
897 IF( tjj.GT.smlnum )
THEN
901 IF( tjj.LT.one )
THEN
902 IF( xj.GT.tjj*bignum )
THEN
907 CALL zdscal( n, rec, x, 1 )
912 x( j ) = zladiv( x( j ), tjjs )
913 ELSE IF( tjj.GT.zero )
THEN
917 IF( xj.GT.tjj*bignum )
THEN
921 rec = ( tjj*bignum ) / xj
922 CALL zdscal( n, rec, x, 1 )
926 x( j ) = zladiv( x( j ), tjjs )
945 x( j ) = zladiv( x( j ), tjjs ) - csumj
947 xmax = max( xmax, cabs1( x( j ) ) )
950 scale = scale / tscal
955 IF( tscal.NE.one )
THEN
956 CALL dscal( n, one / tscal, cnorm, 1 )
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
subroutine ztrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRSV
subroutine zlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
ZLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
subroutine dscal(N, DA, DX, INCX)
DSCAL