262 SUBROUTINE zlasyf_rk( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
272 INTEGER INFO, KB, LDA, LDW, N, NB
276 COMPLEX*16 A( lda, * ), E( * ), W( ldw, * )
282 DOUBLE PRECISION ZERO, ONE
283 parameter ( zero = 0.0d+0, one = 1.0d+0 )
284 DOUBLE PRECISION EIGHT, SEVTEN
285 parameter ( eight = 8.0d+0, sevten = 17.0d+0 )
286 COMPLEX*16 CONE, CZERO
287 parameter ( cone = ( 1.0d+0, 0.0d+0 ),
288 $ czero = ( 0.0d+0, 0.0d+0 ) )
292 INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, K, KK, KW, KKW,
294 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX, SFMIN, DTEMP
295 COMPLEX*16 D11, D12, D21, D22, R1, T, Z
300 DOUBLE PRECISION DLAMCH
301 EXTERNAL lsame, izamax, dlamch
307 INTRINSIC abs, dble, dimag, max, min, sqrt
310 DOUBLE PRECISION CABS1
313 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
321 alpha = ( one+sqrt( sevten ) ) / eight
325 sfmin = dlamch(
'S' )
327 IF( lsame( uplo,
'U' ) )
THEN
349 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
357 CALL zcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
359 $
CALL zgemv(
'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 = izamax( k-1, w( 1, kw ), 1 )
373 colmax = cabs1( w( imax, kw ) )
378 IF( max( absakk, colmax ).EQ.zero )
THEN
385 CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
401 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
420 CALL zcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
421 CALL zcopy( k-imax, a( imax, imax+1 ), lda,
422 $ w( imax+1, kw-1 ), 1 )
425 $
CALL zgemv(
'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 + izamax( k-imax, w( imax+1, kw-1 ),
436 rowmax = cabs1( w( jmax, kw-1 ) )
442 itemp = izamax( imax-1, w( 1, kw-1 ), 1 )
443 dtemp = cabs1( w( itemp, kw-1 ) )
444 IF( dtemp.GT.rowmax )
THEN
454 IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
464 CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
471 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
490 CALL zcopy( 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 zcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
513 CALL zcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
518 CALL zswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
519 CALL zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
528 a( kp, k ) = a( kk, k )
529 CALL zcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
531 CALL zcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
536 CALL zswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
537 CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
541 IF( kstep.EQ.1 )
THEN
551 CALL zcopy( 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 zscal( 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 zgemv(
'No transpose', jj-j+1, n-k, -cone,
639 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
646 $
CALL zgemm(
'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 zcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
682 $
CALL zgemv(
'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 + izamax( n-k, w( k+1, k ), 1 )
696 colmax = cabs1( w( imax, k ) )
701 IF( max( absakk, colmax ).EQ.zero )
THEN
708 CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
724 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
743 CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
744 CALL zcopy( n-imax+1, a( imax, imax ), 1,
745 $ w( imax, k+1 ), 1 )
747 $
CALL zgemv(
'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 + izamax( imax-k, w( k, k+1 ), 1 )
757 rowmax = cabs1( w( jmax, k+1 ) )
763 itemp = imax + izamax( n-imax, w( imax+1, k+1 ), 1)
764 dtemp = cabs1( w( itemp, k+1 ) )
765 IF( dtemp.GT.rowmax )
THEN
775 IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
785 CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
792 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
811 CALL zcopy( 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 zcopy( p-k, a( k, k ), 1, a( p, k ), lda )
830 CALL zcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
835 CALL zswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
836 CALL zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
845 a( kp, k ) = a( kk, k )
846 CALL zcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
847 CALL zcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
851 CALL zswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
852 CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
855 IF( kstep.EQ.1 )
THEN
865 CALL zcopy( 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 zscal( 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 zgemv(
'No transpose', j+jb-jj, k-1, -cone,
952 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
959 $
CALL zgemm(
'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 zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zlasyf_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
ZLASYF_RK computes a partial factorization of a complex symmetric indefinite matrix using bounded Bun...
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL