239 SUBROUTINE zlatrs( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
248 CHARACTER DIAG, NORMIN, TRANS, UPLO
250 DOUBLE PRECISION SCALE
253 DOUBLE PRECISION CNORM( * )
254 COMPLEX*16 A( lda, * ), X( * )
260 DOUBLE PRECISION ZERO, HALF, ONE, TWO
261 parameter ( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0,
265 LOGICAL NOTRAN, NOUNIT, UPPER
266 INTEGER I, IMAX, J, JFIRST, JINC, JLAST
267 DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
269 COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
273 INTEGER IDAMAX, IZAMAX
274 DOUBLE PRECISION DLAMCH, DZASUM
275 COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
276 EXTERNAL lsame, idamax, izamax, dlamch, dzasum, zdotc,
283 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
286 DOUBLE PRECISION CABS1, CABS2
289 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
290 cabs2( zdum ) = abs( dble( zdum ) / 2.d0 ) +
291 $ abs( dimag( zdum ) / 2.d0 )
296 upper = lsame( uplo,
'U' )
297 notran = lsame( trans,
'N' )
298 nounit = lsame( diag,
'N' )
302 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
304 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) .AND. .NOT.
305 $ lsame( trans,
'C' ) )
THEN
307 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag,
'U' ) )
THEN
309 ELSE IF( .NOT.lsame( normin,
'Y' ) .AND. .NOT.
310 $ lsame( normin,
'N' ) )
THEN
312 ELSE IF( n.LT.0 )
THEN
314 ELSE IF( lda.LT.max( 1, n ) )
THEN
318 CALL xerbla(
'ZLATRS', -info )
329 smlnum = dlamch(
'Safe minimum' )
330 bignum = one / smlnum
331 CALL dlabad( smlnum, bignum )
332 smlnum = smlnum / dlamch(
'Precision' )
333 bignum = one / smlnum
336 IF( lsame( normin,
'N' ) )
THEN
345 cnorm( j ) = dzasum( j-1, a( 1, j ), 1 )
352 cnorm( j ) = dzasum( n-j, a( j+1, j ), 1 )
361 imax = idamax( n, cnorm, 1 )
363 IF( tmax.LE.bignum*half )
THEN
366 tscal = half / ( smlnum*tmax )
367 CALL dscal( n, tscal, cnorm, 1 )
375 xmax = max( xmax, cabs2( x( j ) ) )
393 IF( tscal.NE.one )
THEN
405 grow = half / max( xbnd, smlnum )
407 DO 40 j = jfirst, jlast, jinc
417 IF( tjj.GE.smlnum )
THEN
421 xbnd = min( xbnd, min( one, tjj )*grow )
429 IF( tjj+cnorm( j ).GE.smlnum )
THEN
433 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
448 grow = min( one, half / max( xbnd, smlnum ) )
449 DO 50 j = jfirst, jlast, jinc
458 grow = grow*( one / ( one+cnorm( j ) ) )
477 IF( tscal.NE.one )
THEN
489 grow = half / max( xbnd, smlnum )
491 DO 70 j = jfirst, jlast, jinc
500 xj = one + cnorm( j )
501 grow = min( grow, xbnd / xj )
506 IF( tjj.GE.smlnum )
THEN
511 $ xbnd = xbnd*( tjj / xj )
519 grow = min( grow, xbnd )
526 grow = min( one, half / max( xbnd, smlnum ) )
527 DO 80 j = jfirst, jlast, jinc
536 xj = one + cnorm( j )
543 IF( ( grow*tscal ).GT.smlnum )
THEN
548 CALL ztrsv( uplo, trans, diag, n, a, lda, x, 1 )
553 IF( xmax.GT.bignum*half )
THEN
558 scale = ( bignum*half ) / xmax
559 CALL zdscal( n, scale, x, 1 )
569 DO 120 j = jfirst, jlast, jinc
575 tjjs = a( j, j )*tscal
582 IF( tjj.GT.smlnum )
THEN
586 IF( tjj.LT.one )
THEN
587 IF( xj.GT.tjj*bignum )
THEN
592 CALL zdscal( n, rec, x, 1 )
597 x( j ) = zladiv( x( j ), tjjs )
599 ELSE IF( tjj.GT.zero )
THEN
603 IF( xj.GT.tjj*bignum )
THEN
608 rec = ( tjj*bignum ) / xj
609 IF( cnorm( j ).GT.one )
THEN
614 rec = rec / cnorm( j )
616 CALL zdscal( n, rec, x, 1 )
620 x( j ) = zladiv( x( j ), tjjs )
642 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
647 CALL zdscal( n, rec, x, 1 )
650 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
654 CALL zdscal( n, half, x, 1 )
664 CALL zaxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
666 i = izamax( j-1, x, 1 )
667 xmax = cabs1( x( i ) )
675 CALL zaxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
677 i = j + izamax( n-j, x( j+1 ), 1 )
678 xmax = cabs1( x( i ) )
683 ELSE IF( lsame( trans,
'T' ) )
THEN
687 DO 170 j = jfirst, jlast, jinc
694 rec = one / max( xmax, one )
695 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
701 tjjs = a( j, j )*tscal
706 IF( tjj.GT.one )
THEN
710 rec = min( one, rec*tjj )
711 uscal = zladiv( uscal, tjjs )
713 IF( rec.LT.one )
THEN
714 CALL zdscal( n, rec, x, 1 )
721 IF( uscal.EQ.dcmplx( one ) )
THEN
727 csumj = zdotu( j-1, a( 1, j ), 1, x, 1 )
728 ELSE IF( j.LT.n )
THEN
729 csumj = zdotu( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
737 csumj = csumj + ( a( i, j )*uscal )*x( i )
739 ELSE IF( j.LT.n )
THEN
741 csumj = csumj + ( a( i, j )*uscal )*x( i )
746 IF( uscal.EQ.dcmplx( tscal ) )
THEN
751 x( j ) = x( j ) - csumj
754 tjjs = a( j, j )*tscal
764 IF( tjj.GT.smlnum )
THEN
768 IF( tjj.LT.one )
THEN
769 IF( xj.GT.tjj*bignum )
THEN
774 CALL zdscal( n, rec, x, 1 )
779 x( j ) = zladiv( x( j ), tjjs )
780 ELSE IF( tjj.GT.zero )
THEN
784 IF( xj.GT.tjj*bignum )
THEN
788 rec = ( tjj*bignum ) / xj
789 CALL zdscal( n, rec, x, 1 )
793 x( j ) = zladiv( x( j ), tjjs )
812 x( j ) = zladiv( x( j ), tjjs ) - csumj
814 xmax = max( xmax, cabs1( x( j ) ) )
821 DO 220 j = jfirst, jlast, jinc
828 rec = one / max( xmax, one )
829 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
835 tjjs = dconjg( a( j, j ) )*tscal
840 IF( tjj.GT.one )
THEN
844 rec = min( one, rec*tjj )
845 uscal = zladiv( uscal, tjjs )
847 IF( rec.LT.one )
THEN
848 CALL zdscal( n, rec, x, 1 )
855 IF( uscal.EQ.dcmplx( one ) )
THEN
861 csumj = zdotc( j-1, a( 1, j ), 1, x, 1 )
862 ELSE IF( j.LT.n )
THEN
863 csumj = zdotc( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
871 csumj = csumj + ( dconjg( a( i, j ) )*uscal )*
874 ELSE IF( j.LT.n )
THEN
876 csumj = csumj + ( dconjg( a( i, j ) )*uscal )*
882 IF( uscal.EQ.dcmplx( tscal ) )
THEN
887 x( j ) = x( j ) - csumj
890 tjjs = dconjg( a( j, j ) )*tscal
900 IF( tjj.GT.smlnum )
THEN
904 IF( tjj.LT.one )
THEN
905 IF( xj.GT.tjj*bignum )
THEN
910 CALL zdscal( n, rec, x, 1 )
915 x( j ) = zladiv( x( j ), tjjs )
916 ELSE IF( tjj.GT.zero )
THEN
920 IF( xj.GT.tjj*bignum )
THEN
924 rec = ( tjj*bignum ) / xj
925 CALL zdscal( n, rec, x, 1 )
929 x( j ) = zladiv( x( j ), tjjs )
948 x( j ) = zladiv( x( j ), tjjs ) - csumj
950 xmax = max( xmax, cabs1( x( j ) ) )
953 scale = scale / tscal
958 IF( tscal.NE.one )
THEN
959 CALL dscal( n, one / tscal, cnorm, 1 )
subroutine ztrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRSV
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
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 zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY