262 SUBROUTINE clasyf_rk( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
272 INTEGER INFO, KB, LDA, LDW, N, NB
276 COMPLEX A( lda, * ), E( * ), W( ldw, * )
283 parameter ( zero = 0.0e+0, one = 1.0e+0 )
285 parameter ( eight = 8.0e+0, sevten = 17.0e+0 )
287 parameter ( cone = ( 1.0e+0, 0.0e+0 ),
288 $ czero = ( 0.0e+0, 0.0e+0 ) )
292 INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, K, KK, KW, KKW,
294 REAL ABSAKK, ALPHA, COLMAX, ROWMAX, SFMIN, STEMP
295 COMPLEX D11, D12, D21, D22, R1, T, Z
301 EXTERNAL lsame, icamax, slamch
307 INTRINSIC abs, aimag, max, min,
REAL, SQRT
313 cabs1( z ) = abs(
REAL( Z ) ) + abs( AIMAG( z ) )
321 alpha = ( one+sqrt( sevten ) ) / eight
325 sfmin = slamch(
'S' )
327 IF( lsame( uplo,
'U' ) )
THEN
349 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
357 CALL ccopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
359 $
CALL cgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ),
360 $ lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
365 absakk = cabs1( w( k, kw ) )
372 imax = icamax( k-1, w( 1, kw ), 1 )
373 colmax = cabs1( w( imax, kw ) )
378 IF( max( absakk, colmax ).EQ.zero )
THEN
385 CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
401 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
420 CALL ccopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
421 CALL ccopy( k-imax, a( imax, imax+1 ), lda,
422 $ w( imax+1, kw-1 ), 1 )
425 $
CALL cgemv(
'No transpose', k, n-k, -cone,
426 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
427 $ cone, w( 1, kw-1 ), 1 )
434 jmax = imax + icamax( k-imax, w( imax+1, kw-1 ),
436 rowmax = cabs1( w( jmax, kw-1 ) )
442 itemp = icamax( imax-1, w( 1, kw-1 ), 1 )
443 stemp = cabs1( w( itemp, kw-1 ) )
444 IF( stemp.GT.rowmax )
THEN
454 IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
464 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
471 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
490 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
496 IF( .NOT. done )
GOTO 12
508 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
512 CALL ccopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
513 CALL ccopy( p, a( 1, k ), 1, a( 1, p ), 1 )
518 CALL cswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
519 CALL cswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
528 a( kp, k ) = a( kk, k )
529 CALL ccopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
531 CALL ccopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
536 CALL cswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
537 CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
541 IF( kstep.EQ.1 )
THEN
551 CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
553 IF( cabs1( a( k, k ) ).GE.sfmin )
THEN
554 r1 = cone / a( k, k )
555 CALL cscal( k-1, r1, a( 1, k ), 1 )
556 ELSE IF( a( k, k ).NE.czero )
THEN
558 a( ii, k ) = a( ii, k ) / a( k, k )
583 d11 = w( k, kw ) / d12
584 d22 = w( k-1, kw-1 ) / d12
585 t = cone / ( d11*d22-cone )
587 a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
589 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
598 a( k-1, k-1 ) = w( k-1, kw-1 )
600 a( k, k ) = w( k, kw )
601 e( k ) = w( k-1, kw )
612 IF( kstep.EQ.1 )
THEN
632 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
633 jb = min( nb, k-j+1 )
637 DO 40 jj = j, j + jb - 1
638 CALL cgemv(
'No transpose', jj-j+1, n-k, -cone,
639 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
646 $
CALL cgemm(
'No transpose',
'Transpose', j-1, jb,
647 $ n-k, -cone, a( 1, k+1 ), lda, w( j, kw+1 ),
648 $ ldw, cone, a( 1, j ), lda )
672 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
680 CALL ccopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
682 $
CALL cgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
683 $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
688 absakk = cabs1( w( k, k ) )
695 imax = k + icamax( n-k, w( k+1, k ), 1 )
696 colmax = cabs1( w( imax, k ) )
701 IF( max( absakk, colmax ).EQ.zero )
THEN
708 CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
724 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
743 CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
744 CALL ccopy( n-imax+1, a( imax, imax ), 1,
745 $ w( imax, k+1 ), 1 )
747 $
CALL cgemv(
'No transpose', n-k+1, k-1, -cone,
748 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
749 $ cone, w( k, k+1 ), 1 )
756 jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
757 rowmax = cabs1( w( jmax, k+1 ) )
763 itemp = imax + icamax( n-imax, w( imax+1, k+1 ), 1)
764 stemp = cabs1( w( itemp, k+1 ) )
765 IF( stemp.GT.rowmax )
THEN
775 IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
785 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
792 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
811 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
817 IF( .NOT. done )
GOTO 72
825 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
829 CALL ccopy( p-k, a( k, k ), 1, a( p, k ), lda )
830 CALL ccopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
835 CALL cswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
836 CALL cswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
845 a( kp, k ) = a( kk, k )
846 CALL ccopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
847 CALL ccopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
851 CALL cswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
852 CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
855 IF( kstep.EQ.1 )
THEN
865 CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
867 IF( cabs1( a( k, k ) ).GE.sfmin )
THEN
868 r1 = cone / a( k, k )
869 CALL cscal( n-k, r1, a( k+1, k ), 1 )
870 ELSE IF( a( k, k ).NE.czero )
THEN
872 a( ii, k ) = a( ii, k ) / a( k, k )
896 d11 = w( k+1, k+1 ) / d21
897 d22 = w( k, k ) / d21
898 t = cone / ( d11*d22-cone )
900 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
902 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
911 a( k, k ) = w( k, k )
913 a( k+1, k+1 ) = w( k+1, k+1 )
925 IF( kstep.EQ.1 )
THEN
946 jb = min( nb, n-j+1 )
950 DO 100 jj = j, j + jb - 1
951 CALL cgemv(
'No transpose', j+jb-jj, k-1, -cone,
952 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
959 $
CALL cgemm(
'No transpose',
'Transpose', n-j-jb+1, jb,
960 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
961 $ ldw, cone, a( j+jb, j ), lda )
subroutine clasyf_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
CLASYF_RK computes a partial factorization of a complex symmetric indefinite matrix using bounded Bun...
subroutine cscal(N, CA, CX, INCX)
CSCAL
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM