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,
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 )