160 SUBROUTINE chptrf( UPLO, N, AP, IPIV, INFO )
180 parameter ( zero = 0.0e+0, one = 1.0e+0 )
182 parameter ( eight = 8.0e+0, sevten = 17.0e+0 )
186 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
188 REAL ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
190 COMPLEX D12, D21, T, WK, WKM1, WKP1, ZDUM
196 EXTERNAL lsame, icamax, slapy2
202 INTRINSIC abs, aimag, cmplx, conjg, max,
REAL, SQRT
208 cabs1( zdum ) = abs(
REAL( ZDUM ) ) + abs( AIMAG( 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(
'CHPTRF', -info )
228 alpha = ( one+sqrt( sevten ) ) / eight
238 kc = ( n-1 )*n / 2 + 1
251 absakk = abs(
REAL( AP( KC+K-1 ) ) )
257 imax = icamax( 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 ) =
REAL( 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 = icamax( 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(
REAL( AP( KPC+IMAX-1 ) ) ).GE.alpha*
328 CALL cswap( kp-1, ap( knc ), 1, ap( kpc ), 1 )
330 DO 30 j = kp + 1, kk - 1
332 t = conjg( ap( knc+j-1 ) )
333 ap( knc+j-1 ) = conjg( ap( kx ) )
336 ap( kx+kk-1 ) = conjg( ap( kx+kk-1 ) )
337 r1 =
REAL( AP( KNC+KK-1 ) )
338 ap( knc+kk-1 ) =
REAL( AP( KPC+KP-1 ) )
340 IF( kstep.EQ.2 )
THEN
341 ap( kc+k-1 ) =
REAL( AP( KC+K-1 ) )
343 ap( kc+k-2 ) = ap( kc+kp-1 )
347 ap( kc+k-1 ) =
REAL( AP( KC+K-1 ) )
349 $ ap( kc-1 ) =
REAL( AP( KC-1 ) )
354 IF( kstep.EQ.1 )
THEN
366 r1 = one /
REAL( AP( KC+K-1 ) )
367 CALL chpr( uplo, k-1, -r1, ap( kc ), 1, ap )
371 CALL csscal( k-1, r1, ap( kc ), 1 )
388 d = slapy2(
REAL( AP( K-1+( K-1 )*K / 2 ) ),
389 $ aimag( ap( k-1+( k-1 )*k / 2 ) ) )
390 d22 =
REAL( AP( K-1+( K-2 )*( K-1 ) / 2 ) ) / D
391 d11 =
REAL( 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 $ conjg( 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 )*conjg( wk ) -
404 $ ap( i+( k-2 )*( k-1 ) / 2 )*conjg( 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 ) = cmplx(
REAL( AP( J+( J-1 )*
$ J / 2 ) )
418 IF( kstep.EQ.1 )
THEN
453 absakk = abs(
REAL( AP( KC ) ) )
459 imax = k + icamax( n-k, ap( kc+1 ), 1 )
460 colmax = cabs1( ap( kc+imax-k ) )
465 IF( max( absakk, colmax ).EQ.zero )
THEN
472 ap( kc ) =
REAL( AP( KC ) )
474 IF( absakk.GE.alpha*colmax )
THEN
486 DO 70 j = k, imax - 1
487 IF( cabs1( ap( kx ) ).GT.rowmax )
THEN
488 rowmax = cabs1( ap( kx ) )
493 kpc = npp - ( n-imax+1 )*( n-imax+2 ) / 2 + 1
495 jmax = imax + icamax( n-imax, ap( kpc+1 ), 1 )
496 rowmax = max( rowmax, cabs1( ap( kpc+jmax-imax ) ) )
499 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
504 ELSE IF( abs(
REAL( AP( KPC ) ) ).GE.alpha*rowmax ) then
522 $ knc = knc + n - k + 1
529 $
CALL cswap( n-kp, ap( knc+kp-kk+1 ), 1, ap( kpc+1 ),
532 DO 80 j = kk + 1, kp - 1
534 t = conjg( ap( knc+j-kk ) )
535 ap( knc+j-kk ) = conjg( ap( kx ) )
538 ap( knc+kp-kk ) = conjg( ap( knc+kp-kk ) )
539 r1 =
REAL( AP( KNC ) )
540 ap( knc ) =
REAL( AP( KPC ) )
542 IF( kstep.EQ.2 )
THEN
543 ap( kc ) =
REAL( AP( KC ) )
545 ap( kc+1 ) = ap( kc+kp-k )
549 ap( kc ) =
REAL( AP( KC ) )
551 $ ap( knc ) =
REAL( AP( KNC ) )
556 IF( kstep.EQ.1 )
THEN
570 r1 = one /
REAL( AP( KC ) )
571 CALL chpr( uplo, n-k, -r1, ap( kc+1 ), 1,
576 CALL csscal( n-k, r1, ap( kc+1 ), 1 )
597 d = slapy2(
REAL( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ),
598 $ aimag( ap( k+1+( k-1 )*( 2*n-k ) / 2 ) ) )
599 d11 =
REAL( AP( K+1+K*( 2*N-K-1 ) / 2 ) ) / D
600 d22 =
REAL( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D
601 tt = one / ( d11*d22-one )
602 d21 = ap( k+1+( k-1 )*( 2*n-k ) / 2 ) / d
606 wk = d*( d11*ap( j+( k-1 )*( 2*n-k ) / 2 )-d21*
607 $ ap( j+k*( 2*n-k-1 ) / 2 ) )
608 wkp1 = d*( d22*ap( j+k*( 2*n-k-1 ) / 2 )-
609 $ conjg( d21 )*ap( j+( k-1 )*( 2*n-k ) / 2 ) )
611 ap( i+( j-1 )*( 2*n-j ) / 2 ) = ap( i+( j-1 )*
612 $ ( 2*n-j ) / 2 ) - ap( i+( k-1 )*( 2*n-k ) /
613 $ 2 )*conjg( wk ) - ap( i+k*( 2*n-k-1 ) / 2 )*
616 ap( j+( k-1 )*( 2*n-k ) / 2 ) = wk
617 ap( j+k*( 2*n-k-1 ) / 2 ) = wkp1
618 ap( j+( j-1 )*( 2*n-j ) / 2 )
619 $ = cmplx(
REAL( AP( J+( J-1 )*( 2*N-J ) / 2 ) ),
628 IF( kstep.EQ.1 )
THEN
649 subroutine chptrf(UPLO, N, AP, IPIV, INFO)
CHPTRF
subroutine chpr(UPLO, N, ALPHA, X, INCX, AP)
CHPR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine csscal(N, SA, CX, INCX)
CSSCAL