194 INTEGER INFO, KB, LDA, LDW, N, NB
198 COMPLEX A( lda, * ), W( ldw, * )
205 parameter( zero = 0.0e+0, one = 1.0e+0 )
207 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
209 parameter( cone = ( 1.0e+0, 0.0e+0 ),
210 $ czero = ( 0.0e+0, 0.0e+0 ) )
214 INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK,
215 $ kw, kkw, kp, kstep, p, ii
216 REAL ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN
217 COMPLEX D11, D12, D21, D22, R1, T, Z
223 EXTERNAL lsame, icamax, slamch
229 INTRINSIC abs, max, min, sqrt, aimag, real
235 cabs1( z ) = abs(
REAL( Z ) ) + abs( AIMAG( z ) )
243 alpha = ( one+sqrt( sevten ) ) / eight
247 sfmin = slamch(
'S' )
249 IF( lsame( uplo,
'U' ) )
THEN
266 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
274 CALL
ccopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
276 $ CALL
cgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ),
277 $ lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
282 absakk = cabs1( w( k, kw ) )
289 imax = icamax( k-1, w( 1, kw ), 1 )
290 colmax = cabs1( w( imax, kw ) )
295 IF( max( absakk, colmax ).EQ.zero )
THEN
302 CALL
ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
312 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
331 CALL
ccopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
332 CALL
ccopy( k-imax, a( imax, imax+1 ), lda,
333 $ w( imax+1, kw-1 ), 1 )
336 $ CALL
cgemv(
'No transpose', k, n-k, -cone,
337 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
338 $ cone, w( 1, kw-1 ), 1 )
345 jmax = imax + icamax( k-imax, w( imax+1, kw-1 ),
347 rowmax = cabs1( w( jmax, kw-1 ) )
353 itemp = icamax( imax-1, w( 1, kw-1 ), 1 )
354 stemp = cabs1( w( itemp, kw-1 ) )
355 IF( stemp.GT.rowmax )
THEN
365 IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
375 CALL
ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
382 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
401 CALL
ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
407 IF( .NOT. done ) goto 12
419 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
423 CALL
ccopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
424 CALL
ccopy( p, a( 1, k ), 1, a( 1, p ), 1 )
429 CALL
cswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
430 CALL
cswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
439 a( kp, k ) = a( kk, k )
440 CALL
ccopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
442 CALL
ccopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
447 CALL
cswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
448 CALL
cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
452 IF( kstep.EQ.1 )
THEN
462 CALL
ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
464 IF( cabs1( a( k, k ) ).GE.sfmin )
THEN
465 r1 = cone / a( k, k )
466 CALL
cscal( k-1, r1, a( 1, k ), 1 )
467 ELSE IF( a( k, k ).NE.czero )
THEN
469 a( ii, k ) = a( ii, k ) / a( k, k )
489 d11 = w( k, kw ) / d12
490 d22 = w( k-1, kw-1 ) / d12
491 t = cone / ( d11*d22-cone )
493 a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
495 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
502 a( k-1, k-1 ) = w( k-1, kw-1 )
503 a( k-1, k ) = w( k-1, kw )
504 a( k, k ) = w( k, kw )
510 IF( kstep.EQ.1 )
THEN
530 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
531 jb = min( nb, k-j+1 )
535 DO 40 jj = j, j + jb - 1
536 CALL
cgemv(
'No transpose', jj-j+1, n-k, -cone,
537 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
544 $ CALL
cgemm(
'No transpose',
'Transpose', j-1, jb,
545 $ n-k, -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
546 $ cone, a( 1, j ), lda )
567 IF( jp2.NE.jj .AND. j.LE.n )
568 $ CALL
cswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
570 IF( jp1.NE.jj .AND. kstep.EQ.2 )
571 $ CALL
cswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
592 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
600 CALL
ccopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
602 $ CALL
cgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
603 $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
608 absakk = cabs1( w( k, k ) )
615 imax = k + icamax( n-k, w( k+1, k ), 1 )
616 colmax = cabs1( w( imax, k ) )
621 IF( max( absakk, colmax ).EQ.zero )
THEN
628 CALL
ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
638 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
657 CALL
ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
658 CALL
ccopy( n-imax+1, a( imax, imax ), 1,
659 $ w( imax, k+1 ), 1 )
661 $ CALL
cgemv(
'No transpose', n-k+1, k-1, -cone,
662 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
663 $ cone, w( k, k+1 ), 1 )
670 jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
671 rowmax = cabs1( w( jmax, k+1 ) )
677 itemp = imax + icamax( n-imax, w( imax+1, k+1 ), 1)
678 stemp = cabs1( w( itemp, k+1 ) )
679 IF( stemp.GT.rowmax )
THEN
689 IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
699 CALL
ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
706 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
725 CALL
ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
731 IF( .NOT. done ) goto 72
739 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
743 CALL
ccopy( p-k, a( k, k ), 1, a( p, k ), lda )
744 CALL
ccopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
749 CALL
cswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
750 CALL
cswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
759 a( kp, k ) = a( kk, k )
760 CALL
ccopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
761 CALL
ccopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
765 CALL
cswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
766 CALL
cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
769 IF( kstep.EQ.1 )
THEN
779 CALL
ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
781 IF( cabs1( a( k, k ) ).GE.sfmin )
THEN
782 r1 = cone / a( k, k )
783 CALL
cscal( n-k, r1, a( k+1, k ), 1 )
784 ELSE IF( a( k, k ).NE.czero )
THEN
786 a( ii, k ) = a( ii, k ) / a( k, k )
805 d11 = w( k+1, k+1 ) / d21
806 d22 = w( k, k ) / d21
807 t = cone / ( d11*d22-cone )
809 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
811 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
818 a( k, k ) = w( k, k )
819 a( k+1, k ) = w( k+1, k )
820 a( k+1, k+1 ) = w( k+1, k+1 )
826 IF( kstep.EQ.1 )
THEN
847 jb = min( nb, n-j+1 )
851 DO 100 jj = j, j + jb - 1
852 CALL
cgemv(
'No transpose', j+jb-jj, k-1, -cone,
853 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
860 $ CALL
cgemm(
'No transpose',
'Transpose', n-j-jb+1, jb,
861 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ), ldw,
862 $ cone, a( j+jb, j ), lda )
883 IF( jp2.NE.jj .AND. j.GE.1 )
884 $ CALL
cswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
886 IF( jp1.NE.jj .AND. kstep.EQ.2 )
887 $ CALL
cswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
subroutine cscal(N, CA, CX, INCX)
CSCAL
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
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 clasyf_rook(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
CLASYF_ROOK computes a partial factorization of a complex symmetric matrix using the bounded Bunch-Ka...