177 SUBROUTINE dlasyf( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
186 INTEGER INFO, KB, LDA, LDW, N, NB
190 DOUBLE PRECISION A( lda, * ), W( ldw, * )
196 DOUBLE PRECISION ZERO, ONE
197 parameter( zero = 0.0d+0, one = 1.0d+0 )
198 DOUBLE PRECISION EIGHT, SEVTEN
199 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
202 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
204 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
210 EXTERNAL lsame, idamax
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
dcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
249 $ CALL
dgemv(
'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 = idamax( 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
dcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
288 CALL
dcopy( k-imax, a( imax, imax+1 ), lda,
289 $ w( imax+1, kw-1 ), 1 )
291 $ CALL
dgemv(
'No transpose', k, n-k, -one, a( 1, k+1 ),
292 $ lda, w( imax, kw+1 ), ldw, one,
298 jmax = imax + idamax( k-imax, w( imax+1, kw-1 ), 1 )
299 rowmax = abs( w( jmax, kw-1 ) )
301 jmax = idamax( 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
dcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
350 a( kp, kp ) = a( kk, kk )
351 CALL
dcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
354 $ CALL
dcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
362 $ CALL
dswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
364 CALL
dswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
368 IF( kstep.EQ.1 )
THEN
383 CALL
dcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
385 CALL
dscal( 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
dgemv(
'No transpose', jj-j+1, n-k, -one,
486 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
492 CALL
dgemm(
'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
dswap( 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
dcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
545 CALL
dgemv(
'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 + idamax( 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
dcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
584 CALL
dcopy( n-imax+1, a( imax, imax ), 1, w( imax, k+1 ),
586 CALL
dgemv(
'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 + idamax( imax-k, w( k, k+1 ), 1 )
593 rowmax = abs( w( jmax, k+1 ) )
595 jmax = imax + idamax( 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
dcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
640 a( kp, kp ) = a( kk, kk )
641 CALL
dcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
644 $ CALL
dcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
652 $ CALL
dswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
653 CALL
dswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
656 IF( kstep.EQ.1 )
THEN
671 CALL
dcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
674 CALL
dscal( 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
dgemv(
'No transpose', j+jb-jj, k-1, -one,
776 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
783 $ CALL
dgemm(
'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
dswap( j, a( jp, 1 ), lda, a( jj, 1 ), lda )
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dlasyf(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
DLASYF computes a partial factorization of a real symmetric matrix using the Bunch-Kaufman diagonal p...
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV