159 SUBROUTINE zsptrf( UPLO, N, AP, IPIV, INFO )
178 DOUBLE PRECISION ZERO, ONE
179 parameter ( zero = 0.0d+0, one = 1.0d+0 )
180 DOUBLE PRECISION EIGHT, SEVTEN
181 parameter ( eight = 8.0d+0, sevten = 17.0d+0 )
183 parameter ( cone = ( 1.0d+0, 0.0d+0 ) )
187 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
189 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX
190 COMPLEX*16 D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, ZDUM
195 EXTERNAL lsame, izamax
201 INTRINSIC abs, dble, dimag, max, sqrt
204 DOUBLE PRECISION CABS1
207 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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(
'ZSPTRF', -info )
227 alpha = ( one+sqrt( sevten ) ) / eight
237 kc = ( n-1 )*n / 2 + 1
250 absakk = cabs1( ap( kc+k-1 ) )
256 imax = izamax( 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 = izamax( 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 zswap( 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 zspr( uplo, k-1, -r1, ap( kc ), 1, ap )
359 CALL zscal( 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 + izamax( 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 + izamax( 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 zswap( 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 zspr( uplo, n-k, -r1, ap( kc+1 ), 1,
553 CALL zscal( 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 zsptrf(UPLO, N, AP, IPIV, INFO)
ZSPTRF
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zspr(UPLO, N, ALPHA, X, INCX, AP)
ZSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix. ...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL