160 SUBROUTINE zhptrf( UPLO, N, AP, IPIV, INFO )
179 DOUBLE PRECISION zero, one
180 parameter( zero = 0.0d+0, one = 1.0d+0 )
181 DOUBLE PRECISION eight, sevten
182 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
186 INTEGER i, imax,
j, jmax, k, kc, kk, knc, kp, kpc,
188 DOUBLE PRECISION absakk, alpha, colmax, d, d11, d22, r1, rowmax,
190 COMPLEX*16 d12, d21, t, wk, wkm1, wkp1, zdum
202 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, sqrt
205 DOUBLE PRECISION cabs1
208 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
215 upper =
lsame( uplo,
'U' )
216 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
218 ELSE IF( n.LT.0 )
THEN
222 CALL
xerbla(
'ZHPTRF', -info )
228 alpha = ( one+sqrt( sevten ) ) / eight
238 kc = ( n-1 )*n / 2 + 1
251 absakk = abs( dble( ap( kc+k-1 ) ) )
257 imax =
izamax( k-1, ap( kc ), 1 )
258 colmax = cabs1( ap( kc+imax-1 ) )
263 IF( max( absakk, colmax ).EQ.zero )
THEN
270 ap( kc+k-1 ) = dble( ap( kc+k-1 ) )
272 IF( absakk.GE.alpha*colmax )
THEN
284 kx = imax*( imax+1 ) / 2 + imax
285 DO 20
j = imax + 1, k
286 IF( cabs1( ap( kx ) ).GT.rowmax )
THEN
287 rowmax = cabs1( ap( kx ) )
292 kpc = ( imax-1 )*imax / 2 + 1
294 jmax =
izamax( imax-1, ap( kpc ), 1 )
295 rowmax = max( rowmax, cabs1( ap( kpc+jmax-1 ) ) )
298 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
303 ELSE IF( abs( dble( ap( kpc+imax-1 ) ) ).GE.alpha*
328 CALL
zswap( kp-1, ap( knc ), 1, ap( kpc ), 1 )
330 DO 30
j = kp + 1, kk - 1
332 t = dconjg( ap( knc+
j-1 ) )
333 ap( knc+
j-1 ) = dconjg( ap( kx ) )
336 ap( kx+kk-1 ) = dconjg( ap( kx+kk-1 ) )
337 r1 = dble( ap( knc+kk-1 ) )
338 ap( knc+kk-1 ) = dble( ap( kpc+kp-1 ) )
340 IF( kstep.EQ.2 )
THEN
341 ap( kc+k-1 ) = dble( ap( kc+k-1 ) )
343 ap( kc+k-2 ) = ap( kc+kp-1 )
347 ap( kc+k-1 ) = dble( ap( kc+k-1 ) )
349 $ ap( kc-1 ) = dble( ap( kc-1 ) )
354 IF( kstep.EQ.1 )
THEN
366 r1 = one / dble( ap( kc+k-1 ) )
367 CALL
zhpr( uplo, k-1, -r1, ap( kc ), 1, ap )
371 CALL
zdscal( k-1, r1, ap( kc ), 1 )
388 d =
dlapy2( dble( ap( k-1+( k-1 )*k / 2 ) ),
389 $ dimag( ap( k-1+( k-1 )*k / 2 ) ) )
390 d22 = dble( ap( k-1+( k-2 )*( k-1 ) / 2 ) ) / d
391 d11 = dble( ap( k+( k-1 )*k / 2 ) ) / d
392 tt = one / ( d11*d22-one )
393 d12 = ap( k-1+( k-1 )*k / 2 ) / d
396 DO 50
j = k - 2, 1, -1
397 wkm1 = d*( d11*ap(
j+( k-2 )*( k-1 ) / 2 )-
398 $ dconjg( d12 )*ap(
j+( k-1 )*k / 2 ) )
399 wk = d*( d22*ap(
j+( k-1 )*k / 2 )-d12*
400 $ ap(
j+( k-2 )*( k-1 ) / 2 ) )
402 ap( i+(
j-1 )*
j / 2 ) = ap( i+(
j-1 )*
j / 2 ) -
403 $ ap( i+( k-1 )*k / 2 )*dconjg( wk ) -
404 $ ap( i+( k-2 )*( k-1 ) / 2 )*dconjg( wkm1 )
406 ap(
j+( k-1 )*k / 2 ) = wk
407 ap(
j+( k-2 )*( k-1 ) / 2 ) = wkm1
408 ap(
j+(
j-1 )*
j / 2 ) = dcmplx( dble( ap(
j+(
j-
409 $ 1 )*
j / 2 ) ), 0.0d+0 )
419 IF( kstep.EQ.1 )
THEN
454 absakk = abs( dble( ap( kc ) ) )
460 imax = k +
izamax( n-k, ap( kc+1 ), 1 )
461 colmax = cabs1( ap( kc+imax-k ) )
466 IF( max( absakk, colmax ).EQ.zero )
THEN
473 ap( kc ) = dble( ap( kc ) )
475 IF( absakk.GE.alpha*colmax )
THEN
487 DO 70
j = k, imax - 1
488 IF( cabs1( ap( kx ) ).GT.rowmax )
THEN
489 rowmax = cabs1( ap( kx ) )
494 kpc = npp - ( n-imax+1 )*( n-imax+2 ) / 2 + 1
496 jmax = imax +
izamax( n-imax, ap( kpc+1 ), 1 )
497 rowmax = max( rowmax, cabs1( ap( kpc+jmax-imax ) ) )
500 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
505 ELSE IF( abs( dble( ap( kpc ) ) ).GE.alpha*rowmax )
THEN
523 $ knc = knc + n - k + 1
530 $ CALL
zswap( n-kp, ap( knc+kp-kk+1 ), 1, ap( kpc+1 ),
533 DO 80
j = kk + 1, kp - 1
535 t = dconjg( ap( knc+
j-kk ) )
536 ap( knc+
j-kk ) = dconjg( ap( kx ) )
539 ap( knc+kp-kk ) = dconjg( ap( knc+kp-kk ) )
540 r1 = dble( ap( knc ) )
541 ap( knc ) = dble( ap( kpc ) )
543 IF( kstep.EQ.2 )
THEN
544 ap( kc ) = dble( ap( kc ) )
546 ap( kc+1 ) = ap( kc+kp-k )
550 ap( kc ) = dble( ap( kc ) )
552 $ ap( knc ) = dble( ap( knc ) )
557 IF( kstep.EQ.1 )
THEN
571 r1 = one / dble( ap( kc ) )
572 CALL
zhpr( uplo, n-k, -r1, ap( kc+1 ), 1,
577 CALL
zdscal( n-k, r1, ap( kc+1 ), 1 )
598 d =
dlapy2( dble( ap( k+1+( k-1 )*( 2*n-k ) / 2 ) ),
599 $ dimag( ap( k+1+( k-1 )*( 2*n-k ) / 2 ) ) )
600 d11 = dble( ap( k+1+k*( 2*n-k-1 ) / 2 ) ) / d
601 d22 = dble( ap( k+( k-1 )*( 2*n-k ) / 2 ) ) / d
602 tt = one / ( d11*d22-one )
603 d21 = ap( k+1+( k-1 )*( 2*n-k ) / 2 ) / d
607 wk = d*( d11*ap(
j+( k-1 )*( 2*n-k ) / 2 )-d21*
608 $ ap(
j+k*( 2*n-k-1 ) / 2 ) )
609 wkp1 = d*( d22*ap(
j+k*( 2*n-k-1 ) / 2 )-
610 $ dconjg( d21 )*ap(
j+( k-1 )*( 2*n-k ) /
613 ap( i+(
j-1 )*( 2*n-
j ) / 2 ) = ap( i+(
j-1 )*
614 $ ( 2*n-
j ) / 2 ) - ap( i+( k-1 )*( 2*n-k ) /
615 $ 2 )*dconjg( wk ) - ap( i+k*( 2*n-k-1 ) / 2 )*
618 ap(
j+( k-1 )*( 2*n-k ) / 2 ) = wk
619 ap(
j+k*( 2*n-k-1 ) / 2 ) = wkp1
620 ap(
j+(
j-1 )*( 2*n-
j ) / 2 )
621 $ = dcmplx( dble( ap(
j+(
j-1 )*( 2*n-
j ) / 2 ) ),
630 IF( kstep.EQ.1 )
THEN
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
integer function izamax(N, ZX, INCX)
IZAMAX
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zhptrf(UPLO, N, AP, IPIV, INFO)
ZHPTRF
logical function lsame(CA, CB)
LSAME
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine zhpr(UPLO, N, ALPHA, X, INCX, AP)
ZHPR