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
265 INTEGER ICAMAX, ISAMAX
267 COMPLEX CDOTC, CDOTU, CLADIV
268 EXTERNAL lsame, icamax, isamax, scasum, slamch, cdotc,
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
560 CALL csscal( n, scale, x, 1 )
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
594 CALL csscal( n, rec, x, 1 )
599 x( j ) = cladiv( x( j ), tjjs )
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 )
618 CALL csscal( n, rec, x, 1 )
622 x( j ) = cladiv( x( j ), tjjs )
644 IF( cnorm( j ).GT.( bignum-xmax )*rec )
THEN
649 CALL csscal( n, rec, x, 1 )
652 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) )
THEN
656 CALL csscal( n, half, x, 1 )
666 CALL caxpy( j-1, -x( j )*tscal, ap( ip-j+1 ), 1, x,
668 i = icamax( j-1, x, 1 )
669 xmax = cabs1( x( i ) )
678 CALL caxpy( n-j, -x( j )*tscal, ap( ip+1 ), 1,
680 i = j + icamax( n-j, x( j+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
720 CALL csscal( n, rec, x, 1 )
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
780 CALL csscal( n, rec, x, 1 )
785 x( j ) = cladiv( x( j ), tjjs )
786 ELSE IF( tjj.GT.zero )
THEN
790 IF( xj.GT.tjj*bignum )
THEN
794 rec = ( tjj*bignum ) / xj
795 CALL csscal( n, rec, x, 1 )
799 x( j ) = cladiv( x( j ), tjjs )
818 x( j ) = cladiv( x( j ), tjjs ) - csumj
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
858 CALL csscal( n, rec, x, 1 )
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
920 CALL csscal( n, rec, x, 1 )
925 x( j ) = cladiv( x( j ), tjjs )
926 ELSE IF( tjj.GT.zero )
THEN
930 IF( xj.GT.tjj*bignum )
THEN
934 rec = ( tjj*bignum ) / xj
935 CALL csscal( n, rec, x, 1 )
939 x( j ) = cladiv( x( j ), tjjs )
958 x( j ) = cladiv( x( j ), tjjs ) - csumj
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 )
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
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine ctpsv(UPLO, TRANS, DIAG, N, AP, X, INCX)
CTPSV
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
subroutine csscal(N, SA, CX, INCX)
CSSCAL