231 SUBROUTINE clatps( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE,
240 CHARACTER diag, normin, trans, uplo
246 COMPLEX ap( * ),
x( * )
252 REAL zero, half, one, two
253 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
257 LOGICAL notran, nounit, upper
258 INTEGER i, imax, ip,
j, jfirst, jinc, jlast, jlen
259 REAL bignum, grow, rec, smlnum, tjj, tmax, tscal,
261 COMPLEX csumj, tjjs, uscal, zdum
275 INTRINSIC abs, aimag, cmplx, conjg, max, min, real
281 cabs1( zdum ) = abs(
REAL( ZDUM ) ) + abs( aimag( zdum ) )
282 cabs2( zdum ) = abs(
REAL( ZDUM ) / 2. ) +
283 $ abs( aimag( zdum ) / 2. )
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(
'CLATPS', -info )
319 smlnum =
slamch(
'Safe minimum' )
320 bignum = one / smlnum
321 CALL
slabad( smlnum, bignum )
322 smlnum = smlnum /
slamch(
'Precision' )
323 bignum = one / smlnum
326 IF(
lsame( normin,
'N' ) )
THEN
336 cnorm(
j ) =
scasum(
j-1, ap( ip ), 1 )
345 cnorm(
j ) =
scasum( n-
j, ap( ip+1 ), 1 )
355 imax =
isamax( n, cnorm, 1 )
357 IF( tmax.LE.bignum*half )
THEN
360 tscal = half / ( smlnum*tmax )
361 CALL
sscal( 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
ctpsv( 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 110
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
caxpy(
j-1, -
x(
j )*tscal, ap( ip-
j+1 ), 1,
x,
669 xmax = cabs1(
x( i ) )
678 CALL
caxpy( 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 150
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 =
cladiv( uscal, tjjs )
719 IF( rec.LT.one )
THEN
727 IF( uscal.EQ.cmplx( one ) )
THEN
733 csumj =
cdotu(
j-1, ap( ip-
j+1 ), 1,
x, 1 )
734 ELSE IF(
j.LT.n )
THEN
735 csumj =
cdotu( 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.cmplx( 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 190
j = jfirst, jlast, jinc
838 rec = one / max( xmax, one )
839 IF( cnorm(
j ).GT.( bignum-xj )*rec )
THEN
845 tjjs = conjg( ap( ip ) )*tscal
850 IF( tjj.GT.one )
THEN
854 rec = min( one, rec*tjj )
855 uscal =
cladiv( uscal, tjjs )
857 IF( rec.LT.one )
THEN
865 IF( uscal.EQ.cmplx( one ) )
THEN
871 csumj =
cdotc(
j-1, ap( ip-
j+1 ), 1,
x, 1 )
872 ELSE IF(
j.LT.n )
THEN
873 csumj =
cdotc( n-
j, ap( ip+1 ), 1,
x(
j+1 ), 1 )
881 csumj = csumj + ( conjg( ap( ip-
j+i ) )*uscal )*
884 ELSE IF(
j.LT.n )
THEN
886 csumj = csumj + ( conjg( ap( ip+i ) )*uscal )*
892 IF( uscal.EQ.cmplx( tscal ) )
THEN
897 x(
j ) =
x(
j ) - csumj
903 tjjs = conjg( 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
sscal( n, one / tscal, cnorm, 1 )
LOGICAL function lsame(CA, CB)
LSAME
INTEGER function isamax(N, SX, INCX)
ISAMAX
COMPLEX function cladiv(X, Y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
REAL function slamch(CMACH)
SLAMCH
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
subroutine clatps(UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM, INFO)
CLATPS solves a triangular system of equations with the matrix held in packed storage.
subroutine xerbla(SRNAME, INFO)
XERBLA
COMPLEX function cdotu(N, CX, INCX, CY, INCY)
CDOTU
subroutine ctpsv(UPLO, TRANS, DIAG, N, AP, X, INCX)
CTPSV
COMPLEX function cdotc(N, CX, INCX, CY, INCY)
CDOTC
REAL function scasum(N, CX, INCX)
SCASUM
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine csscal(N, SA, CX, INCX)
CSSCAL
INTEGER function icamax(N, CX, INCX)
ICAMAX
subroutine sscal(N, SA, SX, INCX)
SSCAL