178 SUBROUTINE zlasyf( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
187 INTEGER info, kb, lda, ldw, n, nb
191 COMPLEX*16 a( lda, * ), w( ldw, * )
197 DOUBLE PRECISION zero, one
198 parameter( zero = 0.0d+0, one = 1.0d+0 )
199 DOUBLE PRECISION eight, sevten
200 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
202 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
205 INTEGER imax,
j, jb, jj, jmax, jp, k, kk, kkw, kp,
207 DOUBLE PRECISION absakk, alpha, colmax, rowmax
208 COMPLEX*16 d11, d21, d22, r1, t, z
219 INTRINSIC abs, dble, dimag, max, min, sqrt
222 DOUBLE PRECISION cabs1
225 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
233 alpha = ( one+sqrt( sevten ) ) / eight
235 IF(
lsame( uplo,
'U' ) )
THEN
251 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
256 CALL
zcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
258 $ CALL
zgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
259 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
266 absakk = cabs1( w( k, kw ) )
272 imax =
izamax( k-1, w( 1, kw ), 1 )
273 colmax = cabs1( w( imax, kw ) )
278 IF( max( absakk, colmax ).EQ.zero )
THEN
286 IF( absakk.GE.alpha*colmax )
THEN
295 CALL
zcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
296 CALL
zcopy( k-imax, a( imax, imax+1 ), lda,
297 $ w( imax+1, kw-1 ), 1 )
299 $ CALL
zgemv(
'No transpose', k, n-k, -cone,
300 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
301 $ cone, w( 1, kw-1 ), 1 )
306 jmax = imax +
izamax( k-imax, w( imax+1, kw-1 ), 1 )
307 rowmax = cabs1( w( jmax, kw-1 ) )
309 jmax =
izamax( imax-1, w( 1, kw-1 ), 1 )
310 rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
313 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
318 ELSE IF( cabs1( w( imax, kw-1 ) ).GE.alpha*rowmax )
THEN
327 CALL
zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
358 a( kp, kp ) = a( kk, kk )
359 CALL
zcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
362 $ CALL
zcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
370 $ CALL
zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
372 CALL
zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
376 IF( kstep.EQ.1 )
THEN
391 CALL
zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
392 r1 = cone / a( k, k )
393 CALL
zscal( k-1, r1, a( 1, k ), 1 )
440 d11 = w( k, kw ) / d21
441 d22 = w( k-1, kw-1 ) / d21
442 t = cone / ( d11*d22-cone )
450 a(
j, k-1 ) = d21*( d11*w(
j, kw-1 )-w(
j, kw ) )
451 a(
j, k ) = d21*( d22*w(
j, kw )-w(
j, kw-1 ) )
457 a( k-1, k-1 ) = w( k-1, kw-1 )
458 a( k-1, k ) = w( k-1, kw )
459 a( k, k ) = w( k, kw )
467 IF( kstep.EQ.1 )
THEN
487 DO 50
j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
488 jb = min( nb, k-
j+1 )
492 DO 40 jj =
j,
j + jb - 1
493 CALL
zgemv(
'No transpose', jj-
j+1, n-k, -cone,
494 $ a(
j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
500 CALL
zgemm(
'No transpose',
'Transpose',
j-1, jb, n-k,
501 $ -cone, a( 1, k+1 ), lda, w(
j, kw+1 ), ldw,
502 $ cone, a( 1,
j ), lda )
525 IF( jp.NE.jj .AND.
j.LE.n )
526 $ CALL
zswap( n-
j+1, a( jp,
j ), lda, a( jj,
j ), lda )
547 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
552 CALL
zcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
553 CALL
zgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ), lda,
554 $ w( k, 1 ), ldw, cone, w( k, k ), 1 )
561 absakk = cabs1( w( k, k ) )
567 imax = k +
izamax( n-k, w( k+1, k ), 1 )
568 colmax = cabs1( w( imax, k ) )
573 IF( max( absakk, colmax ).EQ.zero )
THEN
581 IF( absakk.GE.alpha*colmax )
THEN
590 CALL
zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
591 CALL
zcopy( n-imax+1, a( imax, imax ), 1, w( imax, k+1 ),
593 CALL
zgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
594 $ lda, w( imax, 1 ), ldw, cone, w( k, k+1 ),
600 jmax = k - 1 +
izamax( imax-k, w( k, k+1 ), 1 )
601 rowmax = cabs1( w( jmax, k+1 ) )
603 jmax = imax +
izamax( n-imax, w( imax+1, k+1 ), 1 )
604 rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
607 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
612 ELSE IF( cabs1( w( imax, k+1 ) ).GE.alpha*rowmax )
THEN
621 CALL
zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
648 a( kp, kp ) = a( kk, kk )
649 CALL
zcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
652 $ CALL
zcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
660 $ CALL
zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
661 CALL
zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
664 IF( kstep.EQ.1 )
THEN
679 CALL
zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
681 r1 = cone / a( k, k )
682 CALL
zscal( n-k, r1, a( k+1, k ), 1 )
730 d11 = w( k+1, k+1 ) / d21
731 d22 = w( k, k ) / d21
732 t = cone / ( d11*d22-cone )
740 a(
j, k ) = d21*( d11*w(
j, k )-w(
j, k+1 ) )
741 a(
j, k+1 ) = d21*( d22*w(
j, k+1 )-w(
j, k ) )
747 a( k, k ) = w( k, k )
748 a( k+1, k ) = w( k+1, k )
749 a( k+1, k+1 ) = w( k+1, k+1 )
757 IF( kstep.EQ.1 )
THEN
778 jb = min( nb, n-
j+1 )
782 DO 100 jj =
j,
j + jb - 1
783 CALL
zgemv(
'No transpose',
j+jb-jj, k-1, -cone,
784 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
791 $ CALL
zgemm(
'No transpose',
'Transpose', n-
j-jb+1, jb,
792 $ k-1, -cone, a(
j+jb, 1 ), lda, w(
j, 1 ),
793 $ ldw, cone, a(
j+jb,
j ), lda )
816 IF( jp.NE.jj .AND.
j.GE.1 )
817 $ CALL
zswap(
j, a( jp, 1 ), lda, a( jj, 1 ), lda )