262 SUBROUTINE slasyf_rk( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
272 INTEGER INFO, KB, LDA, LDW, N, NB
276 REAL 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 )
289 INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, K, KK, KW, KKW,
291 REAL ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
292 $ stemp, r1, rowmax, t, sfmin
298 EXTERNAL lsame, isamax, slamch
304 INTRINSIC abs, max, min, sqrt
312 alpha = ( one+sqrt( sevten ) ) / eight
316 sfmin = slamch(
'S' )
318 IF( lsame( uplo,
'U' ) )
THEN
340 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
348 CALL scopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
350 $
CALL sgemv(
'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 = isamax( k-1, w( 1, kw ), 1 )
364 colmax = abs( w( imax, kw ) )
369 IF( max( absakk, colmax ).EQ.zero )
THEN
376 CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
392 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
411 CALL scopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
412 CALL scopy( k-imax, a( imax, imax+1 ), lda,
413 $ w( imax+1, kw-1 ), 1 )
416 $
CALL sgemv(
'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 + isamax( k-imax, w( imax+1, kw-1 ),
427 rowmax = abs( w( jmax, kw-1 ) )
433 itemp = isamax( imax-1, w( 1, kw-1 ), 1 )
434 stemp = abs( w( itemp, kw-1 ) )
435 IF( stemp.GT.rowmax )
THEN
445 IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
455 CALL scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
462 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
481 CALL scopy( 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 scopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
504 CALL scopy( p, a( 1, k ), 1, a( 1, p ), 1 )
509 CALL sswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
510 CALL sswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
519 a( kp, k ) = a( kk, k )
520 CALL scopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
522 CALL scopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
527 CALL sswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
528 CALL sswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
532 IF( kstep.EQ.1 )
THEN
542 CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
544 IF( abs( a( k, k ) ).GE.sfmin )
THEN
546 CALL sscal( 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 sgemv(
'No transpose', jj-j+1, n-k, -one,
630 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
637 $
CALL sgemm(
'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 scopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
673 $
CALL sgemv(
'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 + isamax( n-k, w( k+1, k ), 1 )
687 colmax = abs( w( imax, k ) )
692 IF( max( absakk, colmax ).EQ.zero )
THEN
699 CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
715 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
734 CALL scopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
735 CALL scopy( n-imax+1, a( imax, imax ), 1,
736 $ w( imax, k+1 ), 1 )
738 $
CALL sgemv(
'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 + isamax( imax-k, w( k, k+1 ), 1 )
748 rowmax = abs( w( jmax, k+1 ) )
754 itemp = imax + isamax( n-imax, w( imax+1, k+1 ), 1)
755 stemp = abs( w( itemp, k+1 ) )
756 IF( stemp.GT.rowmax )
THEN
766 IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
776 CALL scopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
783 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
802 CALL scopy( 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 scopy( p-k, a( k, k ), 1, a( p, k ), lda )
821 CALL scopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
826 CALL sswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
827 CALL sswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
836 a( kp, k ) = a( kk, k )
837 CALL scopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
838 CALL scopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
842 CALL sswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
843 CALL sswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
846 IF( kstep.EQ.1 )
THEN
856 CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
858 IF( abs( a( k, k ) ).GE.sfmin )
THEN
860 CALL sscal( 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 sgemv(
'No transpose', j+jb-jj, k-1, -one,
943 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
950 $
CALL sgemm(
'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 sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine slasyf_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
SLASYF_RK computes a partial factorization of a real symmetric indefinite matrix using bounded Bunch-...
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY