194 INTEGER info, kb, lda, ldw, n, nb
198 DOUBLE PRECISION a( lda, * ), w( ldw, * )
204 DOUBLE PRECISION zero, one
205 parameter( zero = 0.0d+0, one = 1.0d+0 )
206 DOUBLE PRECISION eight, sevten
207 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
211 INTEGER imax, itemp,
j, jb, jj, jmax, jp1, jp2, k, kk,
212 $ kw, kkw, kp, kstep, p, ii
214 DOUBLE PRECISION absakk, alpha, colmax, d11, d12, d21, d22,
215 $ dtemp, r1, rowmax, t, sfmin
227 INTRINSIC abs, max, min, sqrt
235 alpha = ( one+sqrt( sevten ) ) / eight
241 IF(
lsame( uplo,
'U' ) )
THEN
258 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
266 CALL
dcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
268 $ CALL
dgemv(
'No transpose', k, n-k, -one, a( 1, k+1 ),
269 $ lda, w( k, kw+1 ), ldw, one, w( 1, kw ), 1 )
274 absakk = abs( w( k, kw ) )
281 imax =
idamax( k-1, w( 1, kw ), 1 )
282 colmax = abs( w( imax, kw ) )
287 IF( max( absakk, colmax ).EQ.zero )
THEN
294 CALL
dcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
304 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
323 CALL
dcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
324 CALL
dcopy( k-imax, a( imax, imax+1 ), lda,
325 $ w( imax+1, kw-1 ), 1 )
328 $ CALL
dgemv(
'No transpose', k, n-k, -one,
329 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
330 $ one, w( 1, kw-1 ), 1 )
337 jmax = imax +
idamax( k-imax, w( imax+1, kw-1 ),
339 rowmax = abs( w( jmax, kw-1 ) )
345 itemp =
idamax( imax-1, w( 1, kw-1 ), 1 )
346 dtemp = abs( w( itemp, kw-1 ) )
347 IF( dtemp.GT.rowmax )
THEN
357 IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
367 CALL
dcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
374 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
393 CALL
dcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
399 IF( .NOT. done ) goto 12
411 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
415 CALL
dcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
416 CALL
dcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
421 CALL
dswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
422 CALL
dswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
431 a( kp, k ) = a( kk, k )
432 CALL
dcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
434 CALL
dcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
439 CALL
dswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
440 CALL
dswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
444 IF( kstep.EQ.1 )
THEN
454 CALL
dcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
456 IF( abs( a( k, k ) ).GE.sfmin )
THEN
458 CALL
dscal( k-1, r1, a( 1, k ), 1 )
459 ELSE IF( a( k, k ).NE.zero )
THEN
461 a( ii, k ) = a( ii, k ) / a( k, k )
481 d11 = w( k, kw ) / d12
482 d22 = w( k-1, kw-1 ) / d12
483 t = one / ( d11*d22-one )
485 a(
j, k-1 ) = t*( (d11*w(
j, kw-1 )-w(
j, kw ) ) /
487 a(
j, k ) = t*( ( d22*w(
j, kw )-w(
j, kw-1 ) ) /
494 a( k-1, k-1 ) = w( k-1, kw-1 )
495 a( k-1, k ) = w( k-1, kw )
496 a( k, k ) = w( k, kw )
502 IF( kstep.EQ.1 )
THEN
522 DO 50
j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
523 jb = min( nb, k-
j+1 )
527 DO 40 jj =
j,
j + jb - 1
528 CALL
dgemv(
'No transpose', jj-
j+1, n-k, -one,
529 $ a(
j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
536 $ CALL
dgemm(
'No transpose',
'Transpose',
j-1, jb,
537 $ n-k, -one, a( 1, k+1 ), lda, w(
j, kw+1 ), ldw,
538 $ one, a( 1,
j ), lda )
559 IF( jp2.NE.jj .AND.
j.LE.n )
560 $ CALL
dswap( n-
j+1, a( jp2,
j ), lda, a( jj,
j ), lda )
562 IF( jp1.NE.jj .AND. kstep.EQ.2 )
563 $ CALL
dswap( n-
j+1, a( jp1,
j ), lda, a( jj,
j ), lda )
584 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
592 CALL
dcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
594 $ CALL
dgemv(
'No transpose', n-k+1, k-1, -one, a( k, 1 ),
595 $ lda, w( k, 1 ), ldw, one, w( k, k ), 1 )
600 absakk = abs( w( k, k ) )
607 imax = k +
idamax( n-k, w( k+1, k ), 1 )
608 colmax = abs( w( imax, k ) )
613 IF( max( absakk, colmax ).EQ.zero )
THEN
620 CALL
dcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
630 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
649 CALL
dcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
650 CALL
dcopy( n-imax+1, a( imax, imax ), 1,
651 $ w( imax, k+1 ), 1 )
653 $ CALL
dgemv(
'No transpose', n-k+1, k-1, -one,
654 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
655 $ one, w( k, k+1 ), 1 )
662 jmax = k - 1 +
idamax( imax-k, w( k, k+1 ), 1 )
663 rowmax = abs( w( jmax, k+1 ) )
669 itemp = imax +
idamax( n-imax, w( imax+1, k+1 ), 1)
670 dtemp = abs( w( itemp, k+1 ) )
671 IF( dtemp.GT.rowmax )
THEN
681 IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
691 CALL
dcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
698 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
717 CALL
dcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
723 IF( .NOT. done ) goto 72
731 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
735 CALL
dcopy( p-k, a( k, k ), 1, a( p, k ), lda )
736 CALL
dcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
741 CALL
dswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
742 CALL
dswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
751 a( kp, k ) = a( kk, k )
752 CALL
dcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
753 CALL
dcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
757 CALL
dswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
758 CALL
dswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
761 IF( kstep.EQ.1 )
THEN
771 CALL
dcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
773 IF( abs( a( k, k ) ).GE.sfmin )
THEN
775 CALL
dscal( n-k, r1, a( k+1, k ), 1 )
776 ELSE IF( a( k, k ).NE.zero )
THEN
778 a( ii, k ) = a( ii, k ) / a( k, k )
797 d11 = w( k+1, k+1 ) / d21
798 d22 = w( k, k ) / d21
799 t = one / ( d11*d22-one )
801 a(
j, k ) = t*( ( d11*w(
j, k )-w(
j, k+1 ) ) /
803 a(
j, k+1 ) = t*( ( d22*w(
j, k+1 )-w(
j, k ) ) /
810 a( k, k ) = w( k, k )
811 a( k+1, k ) = w( k+1, k )
812 a( k+1, k+1 ) = w( k+1, k+1 )
818 IF( kstep.EQ.1 )
THEN
839 jb = min( nb, n-
j+1 )
843 DO 100 jj =
j,
j + jb - 1
844 CALL
dgemv(
'No transpose',
j+jb-jj, k-1, -one,
845 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
852 $ CALL
dgemm(
'No transpose',
'Transpose', n-
j-jb+1, jb,
853 $ k-1, -one, a(
j+jb, 1 ), lda, w(
j, 1 ), ldw,
854 $ one, a(
j+jb,
j ), lda )
875 IF( jp2.NE.jj .AND.
j.GE.1 )
876 $ CALL
dswap(
j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
878 IF( jp1.NE.jj .AND. kstep.EQ.2 )
879 $ CALL
dswap(
j, a( jp1, 1 ), lda, a( jj, 1 ), lda )