262 SUBROUTINE dlasyf_rk( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
272 INTEGER INFO, KB, LDA, LDW, N, NB
276 DOUBLE PRECISION 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 )
289 INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, K, KK, KW, KKW,
291 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
292 $ dtemp, r1, rowmax, t, sfmin
297 DOUBLE PRECISION DLAMCH
298 EXTERNAL lsame, idamax, dlamch
304 INTRINSIC abs, max, min, sqrt
312 alpha = ( one+sqrt( sevten ) ) / eight
316 sfmin = dlamch(
'S' )
318 IF( lsame( uplo,
'U' ) )
THEN
340 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
348 CALL dcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
350 $
CALL dgemv(
'No transpose', k, n-k, -one, a( 1, k+1 ),
351 $ lda, w( k, kw+1 ), ldw, one, w( 1, kw ), 1 )
356 absakk = abs( w( k, kw ) )
363 imax = idamax( k-1, w( 1, kw ), 1 )
364 colmax = abs( w( imax, kw ) )
369 IF( max( absakk, colmax ).EQ.zero )
THEN
376 CALL dcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
392 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
411 CALL dcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
412 CALL dcopy( k-imax, a( imax, imax+1 ), lda,
413 $ w( imax+1, kw-1 ), 1 )
416 $
CALL dgemv(
'No transpose', k, n-k, -one,
417 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
418 $ one, w( 1, kw-1 ), 1 )
425 jmax = imax + idamax( k-imax, w( imax+1, kw-1 ),
427 rowmax = abs( w( jmax, kw-1 ) )
433 itemp = idamax( imax-1, w( 1, kw-1 ), 1 )
434 dtemp = abs( w( itemp, kw-1 ) )
435 IF( dtemp.GT.rowmax )
THEN
445 IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
455 CALL dcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
462 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
481 CALL dcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
487 IF( .NOT. done )
GOTO 12
499 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
503 CALL dcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
504 CALL dcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
509 CALL dswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
510 CALL dswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
519 a( kp, k ) = a( kk, k )
520 CALL dcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
522 CALL dcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
527 CALL dswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
528 CALL dswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
532 IF( kstep.EQ.1 )
THEN
542 CALL dcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
544 IF( abs( a( k, k ) ).GE.sfmin )
THEN
546 CALL dscal( k-1, r1, a( 1, k ), 1 )
547 ELSE IF( a( k, k ).NE.zero )
THEN
549 a( ii, k ) = a( ii, k ) / a( k, k )
574 d11 = w( k, kw ) / d12
575 d22 = w( k-1, kw-1 ) / d12
576 t = one / ( d11*d22-one )
578 a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
580 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
589 a( k-1, k-1 ) = w( k-1, kw-1 )
591 a( k, k ) = w( k, kw )
592 e( k ) = w( k-1, kw )
603 IF( kstep.EQ.1 )
THEN
623 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
624 jb = min( nb, k-j+1 )
628 DO 40 jj = j, j + jb - 1
629 CALL dgemv(
'No transpose', jj-j+1, n-k, -one,
630 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
637 $
CALL dgemm(
'No transpose',
'Transpose', j-1, jb,
638 $ n-k, -one, a( 1, k+1 ), lda, w( j, kw+1 ),
639 $ ldw, one, a( 1, j ), lda )
663 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
671 CALL dcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
673 $
CALL dgemv(
'No transpose', n-k+1, k-1, -one, a( k, 1 ),
674 $ lda, w( k, 1 ), ldw, one, w( k, k ), 1 )
679 absakk = abs( w( k, k ) )
686 imax = k + idamax( n-k, w( k+1, k ), 1 )
687 colmax = abs( w( imax, k ) )
692 IF( max( absakk, colmax ).EQ.zero )
THEN
699 CALL dcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
715 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
734 CALL dcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
735 CALL dcopy( n-imax+1, a( imax, imax ), 1,
736 $ w( imax, k+1 ), 1 )
738 $
CALL dgemv(
'No transpose', n-k+1, k-1, -one,
739 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
740 $ one, w( k, k+1 ), 1 )
747 jmax = k - 1 + idamax( imax-k, w( k, k+1 ), 1 )
748 rowmax = abs( w( jmax, k+1 ) )
754 itemp = imax + idamax( n-imax, w( imax+1, k+1 ), 1)
755 dtemp = abs( w( itemp, k+1 ) )
756 IF( dtemp.GT.rowmax )
THEN
766 IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
776 CALL dcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
783 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
802 CALL dcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
808 IF( .NOT. done )
GOTO 72
816 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
820 CALL dcopy( p-k, a( k, k ), 1, a( p, k ), lda )
821 CALL dcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
826 CALL dswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
827 CALL dswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
836 a( kp, k ) = a( kk, k )
837 CALL dcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
838 CALL dcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
842 CALL dswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
843 CALL dswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
846 IF( kstep.EQ.1 )
THEN
856 CALL dcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
858 IF( abs( a( k, k ) ).GE.sfmin )
THEN
860 CALL dscal( n-k, r1, a( k+1, k ), 1 )
861 ELSE IF( a( k, k ).NE.zero )
THEN
863 a( ii, k ) = a( ii, k ) / a( k, k )
887 d11 = w( k+1, k+1 ) / d21
888 d22 = w( k, k ) / d21
889 t = one / ( d11*d22-one )
891 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
893 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
902 a( k, k ) = w( k, k )
904 a( k+1, k+1 ) = w( k+1, k+1 )
916 IF( kstep.EQ.1 )
THEN
937 jb = min( nb, n-j+1 )
941 DO 100 jj = j, j + jb - 1
942 CALL dgemv(
'No transpose', j+jb-jj, k-1, -one,
943 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
950 $
CALL dgemm(
'No transpose',
'Transpose', n-j-jb+1, jb,
951 $ k-1, -one, a( j+jb, 1 ), lda, w( j, 1 ),
952 $ ldw, one, a( j+jb, j ), lda )
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dlasyf_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
DLASYF_RK computes a partial factorization of a real symmetric indefinite matrix using bounded Bunch-...
subroutine dscal(N, DA, DX, INCX)
DSCAL