159 SUBROUTINE csptrf( UPLO, N, AP, IPIV, INFO )
179 parameter ( zero = 0.0e+0, one = 1.0e+0 )
181 parameter ( eight = 8.0e+0, sevten = 17.0e+0 )
183 parameter ( cone = ( 1.0e+0, 0.0e+0 ) )
187 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
189 REAL ABSAKK, ALPHA, COLMAX, ROWMAX
190 COMPLEX D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, ZDUM
195 EXTERNAL lsame, icamax
201 INTRINSIC abs, aimag, max,
REAL, SQRT
207 cabs1( zdum ) = abs(
REAL( ZDUM ) ) + abs( AIMAG( zdum ) )
214 upper = lsame( uplo,
'U' )
215 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
217 ELSE IF( n.LT.0 )
THEN
221 CALL xerbla(
'CSPTRF', -info )
227 alpha = ( one+sqrt( sevten ) ) / eight
237 kc = ( n-1 )*n / 2 + 1
250 absakk = cabs1( ap( kc+k-1 ) )
256 imax = icamax( k-1, ap( kc ), 1 )
257 colmax = cabs1( ap( kc+imax-1 ) )
262 IF( max( absakk, colmax ).EQ.zero )
THEN
270 IF( absakk.GE.alpha*colmax )
THEN
279 kx = imax*( imax+1 ) / 2 + imax
280 DO 20 j = imax + 1, k
281 IF( cabs1( ap( kx ) ).GT.rowmax )
THEN
282 rowmax = cabs1( ap( kx ) )
287 kpc = ( imax-1 )*imax / 2 + 1
289 jmax = icamax( imax-1, ap( kpc ), 1 )
290 rowmax = max( rowmax, cabs1( ap( kpc+jmax-1 ) ) )
293 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
298 ELSE IF( cabs1( ap( kpc+imax-1 ) ).GE.alpha*rowmax )
THEN
322 CALL cswap( kp-1, ap( knc ), 1, ap( kpc ), 1 )
324 DO 30 j = kp + 1, kk - 1
327 ap( knc+j-1 ) = ap( kx )
331 ap( knc+kk-1 ) = ap( kpc+kp-1 )
333 IF( kstep.EQ.2 )
THEN
335 ap( kc+k-2 ) = ap( kc+kp-1 )
342 IF( kstep.EQ.1 )
THEN
354 r1 = cone / ap( kc+k-1 )
355 CALL cspr( uplo, k-1, -r1, ap( kc ), 1, ap )
359 CALL cscal( k-1, r1, ap( kc ), 1 )
376 d12 = ap( k-1+( k-1 )*k / 2 )
377 d22 = ap( k-1+( k-2 )*( k-1 ) / 2 ) / d12
378 d11 = ap( k+( k-1 )*k / 2 ) / d12
379 t = cone / ( d11*d22-cone )
382 DO 50 j = k - 2, 1, -1
383 wkm1 = d12*( d11*ap( j+( k-2 )*( k-1 ) / 2 )-
384 $ ap( j+( k-1 )*k / 2 ) )
385 wk = d12*( d22*ap( j+( k-1 )*k / 2 )-
386 $ ap( j+( k-2 )*( k-1 ) / 2 ) )
388 ap( i+( j-1 )*j / 2 ) = ap( i+( j-1 )*j / 2 ) -
389 $ ap( i+( k-1 )*k / 2 )*wk -
390 $ ap( i+( k-2 )*( k-1 ) / 2 )*wkm1
392 ap( j+( k-1 )*k / 2 ) = wk
393 ap( j+( k-2 )*( k-1 ) / 2 ) = wkm1
402 IF( kstep.EQ.1 )
THEN
437 absakk = cabs1( ap( kc ) )
443 imax = k + icamax( n-k, ap( kc+1 ), 1 )
444 colmax = cabs1( ap( kc+imax-k ) )
449 IF( max( absakk, colmax ).EQ.zero )
THEN
457 IF( absakk.GE.alpha*colmax )
THEN
469 DO 70 j = k, imax - 1
470 IF( cabs1( ap( kx ) ).GT.rowmax )
THEN
471 rowmax = cabs1( ap( kx ) )
476 kpc = npp - ( n-imax+1 )*( n-imax+2 ) / 2 + 1
478 jmax = imax + icamax( n-imax, ap( kpc+1 ), 1 )
479 rowmax = max( rowmax, cabs1( ap( kpc+jmax-imax ) ) )
482 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
487 ELSE IF( cabs1( ap( kpc ) ).GE.alpha*rowmax )
THEN
505 $ knc = knc + n - k + 1
512 $
CALL cswap( n-kp, ap( knc+kp-kk+1 ), 1, ap( kpc+1 ),
515 DO 80 j = kk + 1, kp - 1
518 ap( knc+j-kk ) = ap( kx )
522 ap( knc ) = ap( kpc )
524 IF( kstep.EQ.2 )
THEN
526 ap( kc+1 ) = ap( kc+kp-k )
533 IF( kstep.EQ.1 )
THEN
548 CALL cspr( uplo, n-k, -r1, ap( kc+1 ), 1,
553 CALL cscal( n-k, r1, ap( kc+1 ), 1 )
574 d21 = ap( k+1+( k-1 )*( 2*n-k ) / 2 )
575 d11 = ap( k+1+k*( 2*n-k-1 ) / 2 ) / d21
576 d22 = ap( k+( k-1 )*( 2*n-k ) / 2 ) / d21
577 t = cone / ( d11*d22-cone )
581 wk = d21*( d11*ap( j+( k-1 )*( 2*n-k ) / 2 )-
582 $ ap( j+k*( 2*n-k-1 ) / 2 ) )
583 wkp1 = d21*( d22*ap( j+k*( 2*n-k-1 ) / 2 )-
584 $ ap( j+( k-1 )*( 2*n-k ) / 2 ) )
586 ap( i+( j-1 )*( 2*n-j ) / 2 ) = ap( i+( j-1 )*
587 $ ( 2*n-j ) / 2 ) - ap( i+( k-1 )*( 2*n-k ) /
588 $ 2 )*wk - ap( i+k*( 2*n-k-1 ) / 2 )*wkp1
590 ap( j+( k-1 )*( 2*n-k ) / 2 ) = wk
591 ap( j+k*( 2*n-k-1 ) / 2 ) = wkp1
599 IF( kstep.EQ.1 )
THEN
subroutine csptrf(UPLO, N, AP, IPIV, INFO)
CSPTRF
subroutine cscal(N, CA, CX, INCX)
CSCAL
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine cspr(UPLO, N, ALPHA, X, INCX, AP)
CSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix. ...