176 SUBROUTINE clasyf( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
184 INTEGER INFO, KB, LDA, LDW, N, NB
188 COMPLEX A( LDA, * ), W( LDW, * )
195 parameter( zero = 0.0e+0, one = 1.0e+0 )
197 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
199 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
202 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
204 REAL ABSAKK, ALPHA, COLMAX, ROWMAX
205 COMPLEX D11, D21, D22, R1, T, Z
210 EXTERNAL lsame, icamax
216 INTRINSIC abs, aimag, max, min, real, sqrt
222 cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
230 alpha = ( one+sqrt( sevten ) ) / eight
232 IF( lsame( uplo,
'U' ) )
THEN
248 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
253 CALL ccopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
255 $
CALL cgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
256 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
263 absakk = cabs1( w( k, kw ) )
270 imax = icamax( k-1, w( 1, kw ), 1 )
271 colmax = cabs1( w( imax, kw ) )
276 IF( max( absakk, colmax ).EQ.zero )
THEN
284 IF( absakk.GE.alpha*colmax )
THEN
293 CALL ccopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
294 CALL ccopy( k-imax, a( imax, imax+1 ), lda,
295 $ w( imax+1, kw-1 ), 1 )
297 $
CALL cgemv(
'No transpose', k, n-k, -cone,
298 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
299 $ cone, w( 1, kw-1 ), 1 )
304 jmax = imax + icamax( k-imax, w( imax+1, kw-1 ), 1 )
305 rowmax = cabs1( w( jmax, kw-1 ) )
307 jmax = icamax( imax-1, w( 1, kw-1 ), 1 )
308 rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
311 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
316 ELSE IF( cabs1( w( imax, kw-1 ) ).GE.alpha*rowmax )
THEN
325 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
356 a( kp, kp ) = a( kk, kk )
357 CALL ccopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
360 $
CALL ccopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
368 $
CALL cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
370 CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
374 IF( kstep.EQ.1 )
THEN
389 CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
390 r1 = cone / a( k, k )
391 CALL cscal( k-1, r1, a( 1, k ), 1 )
438 d11 = w( k, kw ) / d21
439 d22 = w( k-1, kw-1 ) / d21
440 t = cone / ( d11*d22-cone )
448 a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
449 a( j, k ) = d21*( d22*w( j, kw )-w( j, kw-1 ) )
455 a( k-1, k-1 ) = w( k-1, kw-1 )
456 a( k-1, k ) = w( k-1, kw )
457 a( k, k ) = w( k, kw )
465 IF( kstep.EQ.1 )
THEN
485 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
486 jb = min( nb, k-j+1 )
490 DO 40 jj = j, j + jb - 1
491 CALL cgemv(
'No transpose', jj-j+1, n-k, -cone,
492 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
498 CALL cgemm(
'No transpose',
'Transpose', j-1, jb, n-k,
499 $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
500 $ cone, a( 1, j ), lda )
523 IF( jp.NE.jj .AND. j.LE.n )
524 $
CALL cswap( n-j+1, a( jp, j ), lda, a( jj, j ), lda )
545 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
550 CALL ccopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
551 CALL cgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ), lda,
552 $ w( k, 1 ), ldw, cone, w( k, k ), 1 )
559 absakk = cabs1( w( k, k ) )
566 imax = k + icamax( n-k, w( k+1, k ), 1 )
567 colmax = cabs1( w( imax, k ) )
572 IF( max( absakk, colmax ).EQ.zero )
THEN
580 IF( absakk.GE.alpha*colmax )
THEN
589 CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
590 CALL ccopy( n-imax+1, a( imax, imax ), 1, w( imax, k+1 ),
592 CALL cgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
593 $ lda, w( imax, 1 ), ldw, cone, w( k, k+1 ),
599 jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
600 rowmax = cabs1( w( jmax, k+1 ) )
602 jmax = imax + icamax( n-imax, w( imax+1, k+1 ), 1 )
603 rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
606 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
611 ELSE IF( cabs1( w( imax, k+1 ) ).GE.alpha*rowmax )
THEN
620 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
647 a( kp, kp ) = a( kk, kk )
648 CALL ccopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
651 $
CALL ccopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
659 $
CALL cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
660 CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
663 IF( kstep.EQ.1 )
THEN
678 CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
680 r1 = cone / a( k, k )
681 CALL cscal( n-k, r1, a( k+1, k ), 1 )
729 d11 = w( k+1, k+1 ) / d21
730 d22 = w( k, k ) / d21
731 t = cone / ( d11*d22-cone )
739 a( j, k ) = d21*( d11*w( j, k )-w( j, k+1 ) )
740 a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
746 a( k, k ) = w( k, k )
747 a( k+1, k ) = w( k+1, k )
748 a( k+1, k+1 ) = w( k+1, k+1 )
756 IF( kstep.EQ.1 )
THEN
777 jb = min( nb, n-j+1 )
781 DO 100 jj = j, j + jb - 1
782 CALL cgemv(
'No transpose', j+jb-jj, k-1, -cone,
783 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
790 $
CALL cgemm(
'No transpose',
'Transpose', n-j-jb+1, jb,
791 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
792 $ ldw, cone, a( j+jb, j ), lda )
815 IF( jp.NE.jj .AND. j.GE.1 )
816 $
CALL cswap( j, a( jp, 1 ), lda, a( jj, 1 ), lda )
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
subroutine clasyf(uplo, n, nb, kb, a, lda, ipiv, w, ldw, info)
CLASYF computes a partial factorization of a complex symmetric matrix using the Bunch-Kaufman diagona...
subroutine cscal(n, ca, cx, incx)
CSCAL
subroutine cswap(n, cx, incx, cy, incy)
CSWAP