231 SUBROUTINE zlatps( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE,
240 CHARACTER diag, normin, trans, uplo
242 DOUBLE PRECISION scale
245 DOUBLE PRECISION cnorm( * )
246 COMPLEX*16 ap( * ),
x( * )
252 DOUBLE PRECISION zero, half, one, two
253 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0,
257 LOGICAL notran, nounit, upper
258 INTEGER i, imax, ip,
j, jfirst, jinc, jlast, jlen
259 DOUBLE PRECISION bignum, grow, rec, smlnum, tjj, tmax, tscal,
261 COMPLEX*16 csumj, tjjs, uscal, zdum
275 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
278 DOUBLE PRECISION cabs1, cabs2
281 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
282 cabs2( zdum ) = abs( dble( zdum ) / 2.d0 ) +
283 $ abs( dimag( zdum ) / 2.d0 )
288 upper =
lsame( uplo,
'U' )
289 notran =
lsame( trans,
'N' )
290 nounit =
lsame( diag,
'N' )
294 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
296 ELSE IF( .NOT.notran .AND. .NOT.
lsame( trans,
'T' ) .AND. .NOT.
297 $
lsame( trans,
'C' ) )
THEN
299 ELSE IF( .NOT.nounit .AND. .NOT.
lsame( diag,
'U' ) )
THEN
301 ELSE IF( .NOT.
lsame( normin,
'Y' ) .AND. .NOT.
302 $
lsame( normin,
'N' ) )
THEN
304 ELSE IF( n.LT.0 )
THEN
308 CALL
xerbla(
'ZLATPS', -info )
319 smlnum =
dlamch(
'Safe minimum' )
320 bignum = one / smlnum
321 CALL
dlabad( smlnum, bignum )
322 smlnum = smlnum /
dlamch(
'Precision' )
323 bignum = one / smlnum
326 IF(
lsame( normin,
'N' ) )
THEN
336 cnorm(
j ) =
dzasum(
j-1, ap( ip ), 1 )
345 cnorm(
j ) =
dzasum( n-
j, ap( ip+1 ), 1 )
355 imax =
idamax( n, cnorm, 1 )
357 IF( tmax.LE.bignum*half )
THEN
360 tscal = half / ( smlnum*tmax )
361 CALL
dscal( n, tscal, cnorm, 1 )
369 xmax = max( xmax, cabs2(
x(
j ) ) )
386 IF( tscal.NE.one )
THEN
398 grow = half / max( xbnd, smlnum )
400 ip = jfirst*( jfirst+1 ) / 2
402 DO 40
j = jfirst, jlast, jinc
412 IF( tjj.GE.smlnum )
THEN
416 xbnd = min( xbnd, min( one, tjj )*grow )
424 IF( tjj+cnorm(
j ).GE.smlnum )
THEN
428 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 ip = jfirst*( jfirst+1 ) / 2
490 DO 70
j = jfirst, jlast, jinc
499 xj = one + cnorm(
j )
500 grow = min( grow, xbnd / xj )
505 IF( tjj.GE.smlnum )
THEN
510 $ xbnd = xbnd*( tjj / xj )
520 grow = min( grow, xbnd )
527 grow = min( one, half / max( xbnd, smlnum ) )
528 DO 80
j = jfirst, jlast, jinc
537 xj = one + cnorm(
j )
544 IF( ( grow*tscal ).GT.smlnum )
THEN
549 CALL
ztpsv( uplo, trans, diag, n, ap,
x, 1 )
554 IF( xmax.GT.bignum*half )
THEN
559 scale = ( bignum*half ) / xmax
570 ip = jfirst*( jfirst+1 ) / 2
571 DO 120
j = jfirst, jlast, jinc
577 tjjs = ap( ip )*tscal
584 IF( tjj.GT.smlnum )
THEN
588 IF( tjj.LT.one )
THEN
589 IF( xj.GT.tjj*bignum )
THEN
601 ELSE IF( tjj.GT.zero )
THEN
605 IF( xj.GT.tjj*bignum )
THEN
610 rec = ( tjj*bignum ) / xj
611 IF( cnorm(
j ).GT.one )
THEN
616 rec = rec / cnorm(
j )
644 IF( cnorm(
j ).GT.( bignum-xmax )*rec )
THEN
652 ELSE IF( xj*cnorm(
j ).GT.( bignum-xmax ) )
THEN
666 CALL
zaxpy(
j-1, -
x(
j )*tscal, ap( ip-
j+1 ), 1,
x,
669 xmax = cabs1(
x( i ) )
678 CALL
zaxpy( n-
j, -
x(
j )*tscal, ap( ip+1 ), 1,
681 xmax = cabs1(
x( i ) )
687 ELSE IF(
lsame( trans,
'T' ) )
THEN
691 ip = jfirst*( jfirst+1 ) / 2
693 DO 170
j = jfirst, jlast, jinc
700 rec = one / max( xmax, one )
701 IF( cnorm(
j ).GT.( bignum-xj )*rec )
THEN
707 tjjs = ap( ip )*tscal
712 IF( tjj.GT.one )
THEN
716 rec = min( one, rec*tjj )
717 uscal =
zladiv( uscal, tjjs )
719 IF( rec.LT.one )
THEN
727 IF( uscal.EQ.dcmplx( one ) )
THEN
733 csumj =
zdotu(
j-1, ap( ip-
j+1 ), 1,
x, 1 )
734 ELSE IF(
j.LT.n )
THEN
735 csumj =
zdotu( n-
j, ap( ip+1 ), 1,
x(
j+1 ), 1 )
743 csumj = csumj + ( ap( ip-
j+i )*uscal )*
x( i )
745 ELSE IF(
j.LT.n )
THEN
747 csumj = csumj + ( ap( ip+i )*uscal )*
x(
j+i )
752 IF( uscal.EQ.dcmplx( tscal ) )
THEN
757 x(
j ) =
x(
j ) - csumj
763 tjjs = ap( ip )*tscal
770 IF( tjj.GT.smlnum )
THEN
774 IF( tjj.LT.one )
THEN
775 IF( xj.GT.tjj*bignum )
THEN
786 ELSE IF( tjj.GT.zero )
THEN
790 IF( xj.GT.tjj*bignum )
THEN
794 rec = ( tjj*bignum ) / xj
820 xmax = max( xmax, cabs1(
x(
j ) ) )
829 ip = jfirst*( jfirst+1 ) / 2
831 DO 220
j = jfirst, jlast, jinc
838 rec = one / max( xmax, one )
839 IF( cnorm(
j ).GT.( bignum-xj )*rec )
THEN
845 tjjs = dconjg( ap( ip ) )*tscal
850 IF( tjj.GT.one )
THEN
854 rec = min( one, rec*tjj )
855 uscal =
zladiv( uscal, tjjs )
857 IF( rec.LT.one )
THEN
865 IF( uscal.EQ.dcmplx( one ) )
THEN
871 csumj =
zdotc(
j-1, ap( ip-
j+1 ), 1,
x, 1 )
872 ELSE IF(
j.LT.n )
THEN
873 csumj =
zdotc( n-
j, ap( ip+1 ), 1,
x(
j+1 ), 1 )
881 csumj = csumj + ( dconjg( ap( ip-
j+i ) )*uscal )
884 ELSE IF(
j.LT.n )
THEN
886 csumj = csumj + ( dconjg( ap( ip+i ) )*uscal )*
892 IF( uscal.EQ.dcmplx( tscal ) )
THEN
897 x(
j ) =
x(
j ) - csumj
903 tjjs = dconjg( ap( ip ) )*tscal
910 IF( tjj.GT.smlnum )
THEN
914 IF( tjj.LT.one )
THEN
915 IF( xj.GT.tjj*bignum )
THEN
926 ELSE IF( tjj.GT.zero )
THEN
930 IF( xj.GT.tjj*bignum )
THEN
934 rec = ( tjj*bignum ) / xj
960 xmax = max( xmax, cabs1(
x(
j ) ) )
965 scale = scale / tscal
970 IF( tscal.NE.one )
THEN
971 CALL
dscal( n, one / tscal, cnorm, 1 )
COMPLEX *16 function zladiv(X, Y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
LOGICAL function lsame(CA, CB)
LSAME
COMPLEX *16 function zdotu(N, ZX, INCX, ZY, INCY)
ZDOTU
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
DOUBLE PRECISION function dzasum(N, ZX, INCX)
DZASUM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine ztpsv(UPLO, TRANS, DIAG, N, AP, X, INCX)
ZTPSV
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine zlatps(UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM, INFO)
ZLATPS solves a triangular system of equations with the matrix held in packed storage.
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
COMPLEX *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
INTEGER function izamax(N, ZX, INCX)
IZAMAX
INTEGER function idamax(N, DX, INCX)
IDAMAX