177 SUBROUTINE slasyf( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
186 INTEGER INFO, KB, LDA, LDW, N, NB
190 REAL A( lda, * ), W( ldw, * )
197 parameter( zero = 0.0e+0, one = 1.0e+0 )
199 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
202 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
204 REAL ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
210 EXTERNAL lsame, isamax
216 INTRINSIC abs, max, min, sqrt
224 alpha = ( one+sqrt( sevten ) ) / eight
226 IF( lsame( uplo,
'U' ) )
THEN
242 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
247 CALL
scopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
249 $ CALL
sgemv(
'No transpose', k, n-k, -one, a( 1, k+1 ), lda,
250 $ w( k, kw+1 ), ldw, one, w( 1, kw ), 1 )
257 absakk = abs( w( k, kw ) )
264 imax = isamax( k-1, w( 1, kw ), 1 )
265 colmax = abs( w( imax, kw ) )
270 IF( max( absakk, colmax ).EQ.zero )
THEN
278 IF( absakk.GE.alpha*colmax )
THEN
287 CALL
scopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
288 CALL
scopy( k-imax, a( imax, imax+1 ), lda,
289 $ w( imax+1, kw-1 ), 1 )
291 $ CALL
sgemv(
'No transpose', k, n-k, -one, a( 1, k+1 ),
292 $ lda, w( imax, kw+1 ), ldw, one,
298 jmax = imax + isamax( k-imax, w( imax+1, kw-1 ), 1 )
299 rowmax = abs( w( jmax, kw-1 ) )
301 jmax = isamax( imax-1, w( 1, kw-1 ), 1 )
302 rowmax = max( rowmax, abs( w( jmax, kw-1 ) ) )
305 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
310 ELSE IF( abs( w( imax, kw-1 ) ).GE.alpha*rowmax )
THEN
319 CALL
scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
350 a( kp, kp ) = a( kk, kk )
351 CALL
scopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
354 $ CALL
scopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
362 $ CALL
sswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
364 CALL
sswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
368 IF( kstep.EQ.1 )
THEN
383 CALL
scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
385 CALL
sscal( k-1, r1, a( 1, k ), 1 )
432 d11 = w( k, kw ) / d21
433 d22 = w( k-1, kw-1 ) / d21
434 t = one / ( d11*d22-one )
442 a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
443 a( j, k ) = d21*( d22*w( j, kw )-w( j, kw-1 ) )
449 a( k-1, k-1 ) = w( k-1, kw-1 )
450 a( k-1, k ) = w( k-1, kw )
451 a( k, k ) = w( k, kw )
459 IF( kstep.EQ.1 )
THEN
479 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
480 jb = min( nb, k-j+1 )
484 DO 40 jj = j, j + jb - 1
485 CALL
sgemv(
'No transpose', jj-j+1, n-k, -one,
486 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
492 CALL
sgemm(
'No transpose',
'Transpose', j-1, jb, n-k, -one,
493 $ a( 1, k+1 ), lda, w( j, kw+1 ), ldw, one,
517 IF( jp.NE.jj .AND. j.LE.n )
518 $ CALL
sswap( n-j+1, a( jp, j ), lda, a( jj, j ), lda )
539 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
544 CALL
scopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
545 CALL
sgemv(
'No transpose', n-k+1, k-1, -one, a( k, 1 ), lda,
546 $ w( k, 1 ), ldw, one, w( k, k ), 1 )
553 absakk = abs( w( k, k ) )
560 imax = k + isamax( n-k, w( k+1, k ), 1 )
561 colmax = abs( w( imax, k ) )
566 IF( max( absakk, colmax ).EQ.zero )
THEN
574 IF( absakk.GE.alpha*colmax )
THEN
583 CALL
scopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
584 CALL
scopy( n-imax+1, a( imax, imax ), 1, w( imax, k+1 ),
586 CALL
sgemv(
'No transpose', n-k+1, k-1, -one, a( k, 1 ),
587 $ lda, w( imax, 1 ), ldw, one, w( k, k+1 ), 1 )
592 jmax = k - 1 + isamax( imax-k, w( k, k+1 ), 1 )
593 rowmax = abs( w( jmax, k+1 ) )
595 jmax = imax + isamax( n-imax, w( imax+1, k+1 ), 1 )
596 rowmax = max( rowmax, abs( w( jmax, k+1 ) ) )
599 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
604 ELSE IF( abs( w( imax, k+1 ) ).GE.alpha*rowmax )
THEN
613 CALL
scopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
640 a( kp, kp ) = a( kk, kk )
641 CALL
scopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
644 $ CALL
scopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
652 $ CALL
sswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
653 CALL
sswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
656 IF( kstep.EQ.1 )
THEN
671 CALL
scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
674 CALL
sscal( n-k, r1, a( k+1, k ), 1 )
722 d11 = w( k+1, k+1 ) / d21
723 d22 = w( k, k ) / d21
724 t = one / ( d11*d22-one )
732 a( j, k ) = d21*( d11*w( j, k )-w( j, k+1 ) )
733 a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
739 a( k, k ) = w( k, k )
740 a( k+1, k ) = w( k+1, k )
741 a( k+1, k+1 ) = w( k+1, k+1 )
749 IF( kstep.EQ.1 )
THEN
770 jb = min( nb, n-j+1 )
774 DO 100 jj = j, j + jb - 1
775 CALL
sgemv(
'No transpose', j+jb-jj, k-1, -one,
776 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
783 $ CALL
sgemm(
'No transpose',
'Transpose', n-j-jb+1, jb,
784 $ k-1, -one, a( j+jb, 1 ), lda, w( j, 1 ), ldw,
785 $ one, a( j+jb, j ), lda )
808 IF( jp.NE.jj .AND. j.GE.1 )
809 $ CALL
sswap( j, a( jp, 1 ), lda, a( jj, 1 ), lda )
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine slasyf(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
SLASYF computes a partial factorization of a real symmetric matrix using the Bunch-Kaufman diagonal p...
subroutine sscal(N, SA, SX, INCX)
SSCAL