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
195 DOUBLE PRECISION DLAPY2
196 EXTERNAL lsame, izamax, dlapy2
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
subroutine zhptrf(UPLO, N, AP, IPIV, INFO)
ZHPTRF
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zhpr(UPLO, N, ALPHA, X, INCX, AP)
ZHPR
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL