229 SUBROUTINE dlatps( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE,
238 CHARACTER diag, normin, trans, uplo
240 DOUBLE PRECISION scale
243 DOUBLE PRECISION ap( * ), cnorm( * ), x( * )
249 DOUBLE PRECISION zero, half, one
250 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
253 LOGICAL notran, nounit, upper
254 INTEGER i, imax, ip,
j, jfirst, jinc, jlast, jlen
255 DOUBLE PRECISION bignum, grow, rec, smlnum, sumj, tjj, tjjs,
256 $ tmax, tscal, uscal, xbnd, xj, xmax
268 INTRINSIC abs, max, min
273 upper =
lsame( uplo,
'U' )
274 notran =
lsame( trans,
'N' )
275 nounit =
lsame( diag,
'N' )
279 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
281 ELSE IF( .NOT.notran .AND. .NOT.
lsame( trans,
'T' ) .AND. .NOT.
282 $
lsame( trans,
'C' ) )
THEN
284 ELSE IF( .NOT.nounit .AND. .NOT.
lsame( diag,
'U' ) )
THEN
286 ELSE IF( .NOT.
lsame( normin,
'Y' ) .AND. .NOT.
287 $
lsame( normin,
'N' ) )
THEN
289 ELSE IF( n.LT.0 )
THEN
293 CALL
xerbla(
'DLATPS', -info )
304 smlnum =
dlamch(
'Safe minimum' ) /
dlamch(
'Precision' )
305 bignum = one / smlnum
308 IF(
lsame( normin,
'N' ) )
THEN
318 cnorm(
j ) =
dasum(
j-1, ap( ip ), 1 )
327 cnorm(
j ) =
dasum( n-
j, ap( ip+1 ), 1 )
337 imax =
idamax( n, cnorm, 1 )
339 IF( tmax.LE.bignum )
THEN
342 tscal = one / ( smlnum*tmax )
343 CALL
dscal( n, tscal, cnorm, 1 )
366 IF( tscal.NE.one )
THEN
378 grow = one / max( xbnd, smlnum )
380 ip = jfirst*( jfirst+1 ) / 2
382 DO 30
j = jfirst, jlast, jinc
391 tjj = abs( ap( ip ) )
392 xbnd = min( xbnd, min( one, tjj )*grow )
393 IF( tjj+cnorm(
j ).GE.smlnum )
THEN
397 grow = grow*( tjj / ( tjj+cnorm(
j ) ) )
414 grow = min( one, one / max( xbnd, smlnum ) )
415 DO 40
j = jfirst, jlast, jinc
424 grow = grow*( one / ( one+cnorm(
j ) ) )
443 IF( tscal.NE.one )
THEN
455 grow = one / max( xbnd, smlnum )
457 ip = jfirst*( jfirst+1 ) / 2
459 DO 60
j = jfirst, jlast, jinc
468 xj = one + cnorm(
j )
469 grow = min( grow, xbnd / xj )
473 tjj = abs( ap( ip ) )
475 $ xbnd = xbnd*( tjj / xj )
479 grow = min( grow, xbnd )
486 grow = min( one, one / max( xbnd, smlnum ) )
487 DO 70
j = jfirst, jlast, jinc
496 xj = one + cnorm(
j )
503 IF( ( grow*tscal ).GT.smlnum )
THEN
508 CALL
dtpsv( uplo, trans, diag, n, ap, x, 1 )
513 IF( xmax.GT.bignum )
THEN
518 scale = bignum / xmax
519 CALL
dscal( n, scale, x, 1 )
527 ip = jfirst*( jfirst+1 ) / 2
528 DO 110
j = jfirst, jlast, jinc
534 tjjs = ap( ip )*tscal
541 IF( tjj.GT.smlnum )
THEN
545 IF( tjj.LT.one )
THEN
546 IF( xj.GT.tjj*bignum )
THEN
551 CALL
dscal( n, rec, x, 1 )
556 x(
j ) = x(
j ) / tjjs
558 ELSE IF( tjj.GT.zero )
THEN
562 IF( xj.GT.tjj*bignum )
THEN
567 rec = ( tjj*bignum ) / xj
568 IF( cnorm(
j ).GT.one )
THEN
573 rec = rec / cnorm(
j )
575 CALL
dscal( n, rec, x, 1 )
579 x(
j ) = x(
j ) / tjjs
601 IF( cnorm(
j ).GT.( bignum-xmax )*rec )
THEN
606 CALL
dscal( n, rec, x, 1 )
609 ELSE IF( xj*cnorm(
j ).GT.( bignum-xmax ) )
THEN
613 CALL
dscal( n, half, x, 1 )
623 CALL
daxpy(
j-1, -x(
j )*tscal, ap( ip-
j+1 ), 1, x,
635 CALL
daxpy( n-
j, -x(
j )*tscal, ap( ip+1 ), 1,
648 ip = jfirst*( jfirst+1 ) / 2
650 DO 160
j = jfirst, jlast, jinc
657 rec = one / max( xmax, one )
658 IF( cnorm(
j ).GT.( bignum-xj )*rec )
THEN
664 tjjs = ap( ip )*tscal
669 IF( tjj.GT.one )
THEN
673 rec = min( one, rec*tjj )
676 IF( rec.LT.one )
THEN
677 CALL
dscal( n, rec, x, 1 )
684 IF( uscal.EQ.one )
THEN
690 sumj =
ddot(
j-1, ap( ip-
j+1 ), 1, x, 1 )
691 ELSE IF(
j.LT.n )
THEN
692 sumj =
ddot( n-
j, ap( ip+1 ), 1, x(
j+1 ), 1 )
700 sumj = sumj + ( ap( ip-
j+i )*uscal )*x( i )
702 ELSE IF(
j.LT.n )
THEN
704 sumj = sumj + ( ap( ip+i )*uscal )*x(
j+i )
709 IF( uscal.EQ.tscal )
THEN
714 x(
j ) = x(
j ) - sumj
720 tjjs = ap( ip )*tscal
727 IF( tjj.GT.smlnum )
THEN
731 IF( tjj.LT.one )
THEN
732 IF( xj.GT.tjj*bignum )
THEN
737 CALL
dscal( n, rec, x, 1 )
742 x(
j ) = x(
j ) / tjjs
743 ELSE IF( tjj.GT.zero )
THEN
747 IF( xj.GT.tjj*bignum )
THEN
751 rec = ( tjj*bignum ) / xj
752 CALL
dscal( n, rec, x, 1 )
756 x(
j ) = x(
j ) / tjjs
775 x(
j ) = x(
j ) / tjjs - sumj
777 xmax = max( xmax, abs( x(
j ) ) )
782 scale = scale / tscal
787 IF( tscal.NE.one )
THEN
788 CALL
dscal( n, one / tscal, cnorm, 1 )