252 COMPLEX*16 A( LDA, * ), E( * )
258 DOUBLE PRECISION ZERO, ONE
259 parameter( zero = 0.0d+0, one = 1.0d+0 )
260 DOUBLE PRECISION EIGHT, SEVTEN
261 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
262 COMPLEX*16 CONE, CZERO
263 parameter( cone = ( 1.0d+0, 0.0d+0 ),
264 $ czero = ( 0.0d+0, 0.0d+0 ) )
268 INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
270 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX, DTEMP, SFMIN
271 COMPLEX*16 D11, D12, D21, D22, T, WK, WKM1, WKP1, Z
276 DOUBLE PRECISION DLAMCH
277 EXTERNAL lsame, izamax, dlamch
283 INTRINSIC abs, max, sqrt, dimag, dble
286 DOUBLE PRECISION CABS1
289 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
296 upper = lsame( uplo,
'U' )
297 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
299 ELSE IF( n.LT.0 )
THEN
301 ELSE IF( lda.LT.max( 1, n ) )
THEN
305 CALL xerbla(
'ZSYTF2_RK', -info )
311 alpha = ( one+sqrt( sevten ) ) / eight
315 sfmin = dlamch(
'S' )
342 absakk = cabs1( a( k, k ) )
349 imax = izamax( k-1, a( 1, k ), 1 )
350 colmax = cabs1( a( imax, k ) )
355 IF( (max( absakk, colmax ).EQ.zero) )
THEN
375 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
396 jmax = imax + izamax( k-imax, a( imax, imax+1 ),
398 rowmax = cabs1( a( imax, jmax ) )
404 itemp = izamax( imax-1, a( 1, imax ), 1 )
405 dtemp = cabs1( a( itemp, imax ) )
406 IF( dtemp.GT.rowmax )
THEN
415 IF( .NOT.( cabs1( a( imax, imax ) ).LT.alpha*rowmax ))
427 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) )
THEN
446 IF( .NOT. done )
GOTO 12
454 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
460 $
CALL zswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
462 $
CALL zswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
465 a( k, k ) = a( p, p )
472 $
CALL zswap( n-k, a( k, k+1 ), lda, a( p, k+1 ), lda )
485 $
CALL zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
486 IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
487 $
CALL zswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
490 a( kk, kk ) = a( kp, kp )
492 IF( kstep.EQ.2 )
THEN
494 a( k-1, k ) = a( kp, k )
502 $
CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
509 IF( kstep.EQ.1 )
THEN
522 IF( cabs1( a( k, k ) ).GE.sfmin )
THEN
528 d11 = cone / a( k, k )
529 CALL zsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
533 CALL zscal( k-1, d11, a( 1, k ), 1 )
540 a( ii, k ) = a( ii, k ) / d11
548 CALL zsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
576 d22 = a( k-1, k-1 ) / d12
577 d11 = a( k, k ) / d12
578 t = cone / ( d11*d22-cone )
580 DO 30 j = k - 2, 1, -1
582 wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
583 wk = t*( d22*a( j, k )-a( j, k-1 ) )
586 a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
587 $ ( a( i, k-1 ) / d12 )*wkm1
593 a( j, k-1 ) = wkm1 / d12
614 IF( kstep.EQ.1 )
THEN
652 absakk = cabs1( a( k, k ) )
659 imax = k + izamax( n-k, a( k+1, k ), 1 )
660 colmax = cabs1( a( imax, k ) )
665 IF( ( max( absakk, colmax ).EQ.zero ) )
THEN
685 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
706 jmax = k - 1 + izamax( imax-k, a( imax, k ), lda )
707 rowmax = cabs1( a( imax, jmax ) )
713 itemp = imax + izamax( n-imax, a( imax+1, imax ),
715 dtemp = cabs1( a( itemp, imax ) )
716 IF( dtemp.GT.rowmax )
THEN
725 IF( .NOT.( cabs1( a( imax, imax ) ).LT.alpha*rowmax ))
737 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) )
THEN
756 IF( .NOT. done )
GOTO 42
764 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
770 $
CALL zswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
772 $
CALL zswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
774 a( k, k ) = a( p, p )
781 $
CALL zswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
794 $
CALL zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
795 IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
796 $
CALL zswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
799 a( kk, kk ) = a( kp, kp )
801 IF( kstep.EQ.2 )
THEN
803 a( k+1, k ) = a( kp, k )
811 $
CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
817 IF( kstep.EQ.1 )
THEN
830 IF( cabs1( a( k, k ) ).GE.sfmin )
THEN
836 d11 = cone / a( k, k )
837 CALL zsyr( uplo, n-k, -d11, a( k+1, k ), 1,
838 $ a( k+1, k+1 ), lda )
842 CALL zscal( n-k, d11, a( k+1, k ), 1 )
849 a( ii, k ) = a( ii, k ) / d11
857 CALL zsyr( uplo, n-k, -d11, a( k+1, k ), 1,
858 $ a( k+1, k+1 ), lda )
887 d11 = a( k+1, k+1 ) / d21
888 d22 = a( k, k ) / d21
889 t = cone / ( d11*d22-cone )
895 wk = t*( d11*a( j, k )-a( j, k+1 ) )
896 wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
901 a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
902 $ ( a( i, k+1 ) / d21 )*wkp1
908 a( j, k+1 ) = wkp1 / d21
929 IF( kstep.EQ.1 )
THEN
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
subroutine zsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
ZSYR performs the symmetric rank-1 update of a complex symmetric matrix.
subroutine zsytf2_rk(UPLO, N, A, LDA, E, IPIV, INFO)
ZSYTF2_RK computes the factorization of a complex symmetric indefinite matrix using the bounded Bunch...