1 SUBROUTINE zlatrs( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
10 CHARACTER DIAG, NORMIN, TRANS, UPLO
12 DOUBLE PRECISION SCALE
15 DOUBLE PRECISION CNORM( * )
16 COMPLEX*16 A( LDA, * ), X( * )
174 DOUBLE PRECISION ZERO, HALF, ONE, TWO
175 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0,
179 LOGICAL NOTRAN, NOUNIT, UPPER
180 INTEGER I, IMAX, J, JFIRST, JINC, JLAST
181 DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
183 COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
187 INTEGER IDAMAX, IZAMAX
188 DOUBLE PRECISION DLAMCH, DZASUM
189 COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
190 EXTERNAL lsame, idamax, izamax, dlamch, dzasum, zdotc,
197 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
200 DOUBLE PRECISION CABS1, CABS2
203 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
204 cabs2( zdum ) = abs( dble( zdum ) / 2.d0 ) +
205 $ abs( dimag( zdum ) / 2.d0 )
210 upper = lsame( uplo,
'U' )
211 notran = lsame( trans,
'N' )
212 nounit = lsame( diag,
'N' )
216 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
218 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) .AND. .NOT.
219 $ lsame( trans,
'C' ) )
THEN
221 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag,
'U' ) )
THEN
223 ELSE IF( .NOT.lsame( normin,
'Y' ) .AND. .NOT.
224 $ lsame( normin,
'N' ) )
THEN
226 ELSE IF( n.LT.0 )
THEN
228 ELSE IF( lda.LT.max( 1, n ) )
THEN
232 CALL xerbla(
'ZLATRS', -info )
243 smlnum = dlamch(
'Safe minimum' )
244 bignum = one / smlnum
245 CALL dlabad( smlnum, bignum )
246 smlnum = smlnum / dlamch(
'Precision' )
247 bignum = one / smlnum
250 IF( lsame( normin,
'N' ) )
THEN
259 cnorm( j ) = dzasum( j-1, a( 1, j ), 1 )
266 cnorm( j ) = dzasum( n-j, a( j+1, j ), 1 )
275 imax = idamax( n, cnorm, 1 )
277 IF( tmax.LE.bignum*half )
THEN
280 tscal = half / ( smlnum*tmax )
281 CALL dscal( n, tscal, cnorm, 1 )
289 xmax = max( xmax, cabs2( x( j ) ) )
307 IF( tscal.NE.one )
THEN
319 grow = half / max( xbnd, smlnum )
321 DO 40 j = jfirst, jlast, jinc
331 IF( tjj.GE.smlnum )
THEN
335 xbnd = min( xbnd, min( one, tjj )*grow )
343 IF( tjj+cnorm( j ).GE.smlnum )
THEN
347 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
362 grow = min( one, half / max( xbnd, smlnum ) )
363 DO 50 j = jfirst, jlast, jinc
372 grow = grow*( one / ( one+cnorm( j ) ) )
391 IF( tscal.NE.one )
THEN
403 grow = half / max( xbnd, smlnum )
405 DO 70 j = jfirst, jlast, jinc
414 xj = one + cnorm( j )
415 grow = min( grow, xbnd / xj )
420 IF( tjj.GE.smlnum )
THEN
425 $ xbnd = xbnd*( tjj / xj )
433 grow = min( grow, xbnd )
440 grow = min( one, half / max( xbnd, smlnum ) )
441 DO 80 j = jfirst, jlast, jinc
450 xj = one + cnorm( j )
457 IF( ( grow*tscal ).GT.smlnum )
THEN
462 CALL ztrsv( uplo, trans, diag, n, a, lda, x, 1 )
467 IF( xmax.GT.bignum*half )
THEN
472 scale = ( bignum*half ) / xmax
473 CALL zdscal( n, scale, x, 1 )
483 DO 120 j = jfirst, jlast, jinc
489 tjjs = a( j, j )*tscal
496 IF( tjj.GT.smlnum )
THEN
500 IF( tjj.LT.one )
THEN
501 IF( xj.GT.tjj*bignum )
THEN
506 CALL zdscal( n, rec, x, 1 )
511 x( j ) = zladiv( x( j ), tjjs )
513 ELSE IF( tjj.GT.zero )
THEN
517 IF( xj.GT.tjj*bignum )
THEN
522 rec = ( tjj*bignum ) / xj
523 IF( cnorm( j ).GT.one )
THEN
528 rec = rec / cnorm( j )
530 CALL zdscal( n, rec, x, 1 )
534 x( j ) = zladiv( x( j ), tjjs )
556 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
561 CALL zdscal( n, rec, x, 1 )
564 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
568 CALL zdscal( n, half, x, 1 )
578 CALL zaxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
580 i = izamax( j-1, x, 1 )
581 xmax = cabs1( x( i ) )
589 CALL zaxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
591 i = j + izamax( n-j, x( j+1 ), 1 )
592 xmax = cabs1( x( i ) )
597 ELSE IF( lsame( trans,
'T' ) )
THEN
601 DO 170 j = jfirst, jlast, jinc
608 rec = one / max( xmax, one )
609 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
615 tjjs = a( j, j )*tscal
620 IF( tjj.GT.one )
THEN
624 rec = min( one, rec*tjj )
625 uscal = zladiv( uscal, tjjs )
627 IF( rec.LT.one )
THEN
628 CALL zdscal( n, rec, x, 1 )
635 IF( uscal.EQ.dcmplx( one ) )
THEN
641 csumj = zdotu( j-1, a( 1, j ), 1, x, 1 )
642 ELSE IF( j.LT.n )
THEN
643 csumj = zdotu( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
651 csumj = csumj + ( a( i, j )*uscal )*x( i )
653 ELSE IF( j.LT.n )
THEN
655 csumj = csumj + ( a( i, j )*uscal )*x( i )
660 IF( uscal.EQ.dcmplx( tscal ) )
THEN
665 x( j ) = x( j ) - csumj
668 tjjs = a( j, j )*tscal
678 IF( tjj.GT.smlnum )
THEN
682 IF( tjj.LT.one )
THEN
683 IF( xj.GT.tjj*bignum )
THEN
688 CALL zdscal( n, rec, x, 1 )
693 x( j ) = zladiv( x( j ), tjjs )
694 ELSE IF( tjj.GT.zero )
THEN
698 IF( xj.GT.tjj*bignum )
THEN
702 rec = ( tjj*bignum ) / xj
703 CALL zdscal( n, rec, x, 1 )
707 x( j ) = zladiv( x( j ), tjjs )
726 x( j ) = zladiv( x( j ), tjjs ) - csumj
728 xmax = max( xmax, cabs1( x( j ) ) )
735 DO 220 j = jfirst, jlast, jinc
742 rec = one / max( xmax, one )
743 IF( cnorm( j ).GT.( bignum-xj )*rec )
THEN
749 tjjs = dconjg( a( j, j ) )*tscal
754 IF( tjj.GT.one )
THEN
758 rec = min( one, rec*tjj )
759 uscal = zladiv( uscal, tjjs )
761 IF( rec.LT.one )
THEN
762 CALL zdscal( n, rec, x, 1 )
769 IF( uscal.EQ.dcmplx( one ) )
THEN
775 csumj = zdotc( j-1, a( 1, j ), 1, x, 1 )
776 ELSE IF( j.LT.n )
THEN
777 csumj = zdotc( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
785 csumj = csumj + ( dconjg( a( i, j ) )*uscal )*
788 ELSE IF( j.LT.n )
THEN
790 csumj = csumj + ( dconjg( a( i, j ) )*uscal )*
796 IF( uscal.EQ.dcmplx( tscal ) )
THEN
801 x( j ) = x( j ) - csumj
804 tjjs = dconjg( a( j, j ) )*tscal
814 IF( tjj.GT.smlnum )
THEN
818 IF( tjj.LT.one )
THEN
819 IF( xj.GT.tjj*bignum )
THEN
824 CALL zdscal( n, rec, x, 1 )
829 x( j ) = zladiv( x( j ), tjjs )
830 ELSE IF( tjj.GT.zero )
THEN
834 IF( xj.GT.tjj*bignum )
THEN
838 rec = ( tjj*bignum ) / xj
839 CALL zdscal( n, rec, x, 1 )
843 x( j ) = zladiv( x( j ), tjjs )
862 x( j ) = zladiv( x( j ), tjjs ) - csumj
864 xmax = max( xmax, cabs1( x( j ) ) )
867 scale = scale / tscal
872 IF( tscal.NE.one )
THEN
873 CALL dscal( n, one / tscal, cnorm, 1 )
subroutine dlabad(SMALL, LARGE)
subroutine dscal(n, da, dx, incx)
subroutine xerbla(SRNAME, INFO)
subroutine zaxpy(n, za, zx, incx, zy, incy)
subroutine zdscal(n, da, zx, incx)
subroutine zlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
subroutine ztrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)