154 SUBROUTINE slasd4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
166 REAL D( * ), DELTA( * ), WORK( * ), Z( * )
173 parameter ( maxit = 400 )
174 REAL ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
175 parameter ( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0,
176 $ three = 3.0e+0, four = 4.0e+0, eight = 8.0e+0,
180 LOGICAL ORGATI, SWTCH, SWTCH3, GEOMAVG
181 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
182 REAL A, B, C, DELSQ, DELSQ2, SQ2, DPHI, DPSI, DTIIM,
183 $ dtiip, dtipsq, dtisq, dtnsq, dtnsq1, dw, eps,
184 $ erretm, eta, phi, prew, psi, rhoinv, sglb,
185 $ sgub, tau, tau2, temp, temp1, temp2, w
188 REAL DD( 3 ), ZZ( 3 )
198 INTRINSIC abs, max, min, sqrt
212 sigma = sqrt( d( 1 )*d( 1 )+rho*z( 1 )*z( 1 ) )
218 CALL slasd5( i, d, z, delta, rho, sigma, work )
224 eps = slamch(
'Epsilon' )
244 temp1 = temp / ( d( n )+sqrt( d( n )*d( n )+temp ) )
246 work( j ) = d( j ) + d( n ) + temp1
247 delta( j ) = ( d( j )-d( n ) ) - temp1
252 psi = psi + z( j )*z( j ) / ( delta( j )*work( j ) )
256 w = c + z( ii )*z( ii ) / ( delta( ii )*work( ii ) ) +
257 $ z( n )*z( n ) / ( delta( n )*work( n ) )
260 temp1 = sqrt( d( n )*d( n )+rho )
261 temp = z( n-1 )*z( n-1 ) / ( ( d( n-1 )+temp1 )*
262 $ ( d( n )-d( n-1 )+rho / ( d( n )+temp1 ) ) ) +
263 $ z( n )*z( n ) / rho
271 delsq = ( d( n )-d( n-1 ) )*( d( n )+d( n-1 ) )
272 a = -c*delsq + z( n-1 )*z( n-1 ) + z( n )*z( n )
273 b = z( n )*z( n )*delsq
275 tau2 = two*b / ( sqrt( a*a+four*b*c )-a )
277 tau2 = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
279 tau = tau2 / ( d( n )+sqrt( d( n )*d( n )+tau2 ) )
286 delsq = ( d( n )-d( n-1 ) )*( d( n )+d( n-1 ) )
287 a = -c*delsq + z( n-1 )*z( n-1 ) + z( n )*z( n )
288 b = z( n )*z( n )*delsq
294 tau2 = two*b / ( sqrt( a*a+four*b*c )-a )
296 tau2 = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
298 tau = tau2 / ( d( n )+sqrt( d( n )*d( n )+tau2 ) )
312 delta( j ) = ( d( j )-d( n ) ) - tau
313 work( j ) = d( j ) + d( n ) + tau
322 temp = z( j ) / ( delta( j )*work( j ) )
323 psi = psi + z( j )*temp
324 dpsi = dpsi + temp*temp
325 erretm = erretm + psi
327 erretm = abs( erretm )
331 temp = z( n ) / ( delta( n )*work( n ) )
334 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv
337 w = rhoinv + phi + psi
341 IF( abs( w ).LE.eps*erretm )
THEN
348 dtnsq1 = work( n-1 )*delta( n-1 )
349 dtnsq = work( n )*delta( n )
350 c = w - dtnsq1*dpsi - dtnsq*dphi
351 a = ( dtnsq+dtnsq1 )*w - dtnsq*dtnsq1*( dpsi+dphi )
356 eta = rho - sigma*sigma
357 ELSE IF( a.GE.zero )
THEN
358 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
360 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
370 $ eta = -w / ( dpsi+dphi )
375 eta = eta / ( sigma+sqrt( eta+sigma*sigma ) )
380 delta( j ) = delta( j ) - eta
381 work( j ) = work( j ) + eta
390 temp = z( j ) / ( work( j )*delta( j ) )
391 psi = psi + z( j )*temp
392 dpsi = dpsi + temp*temp
393 erretm = erretm + psi
395 erretm = abs( erretm )
399 tau2 = work( n )*delta( n )
403 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv
406 w = rhoinv + phi + psi
412 DO 90 niter = iter, maxit
416 IF( abs( w ).LE.eps*erretm )
THEN
422 dtnsq1 = work( n-1 )*delta( n-1 )
423 dtnsq = work( n )*delta( n )
424 c = w - dtnsq1*dpsi - dtnsq*dphi
425 a = ( dtnsq+dtnsq1 )*w - dtnsq1*dtnsq*( dpsi+dphi )
428 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
430 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
440 $ eta = -w / ( dpsi+dphi )
445 eta = eta / ( sigma+sqrt( eta+sigma*sigma ) )
450 delta( j ) = delta( j ) - eta
451 work( j ) = work( j ) + eta
460 temp = z( j ) / ( work( j )*delta( j ) )
461 psi = psi + z( j )*temp
462 dpsi = dpsi + temp*temp
463 erretm = erretm + psi
465 erretm = abs( erretm )
469 tau2 = work( n )*delta( n )
473 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv
476 w = rhoinv + phi + psi
495 delsq = ( d( ip1 )-d( i ) )*( d( ip1 )+d( i ) )
497 sq2=sqrt( ( d( i )*d( i )+d( ip1 )*d( ip1 ) ) / two )
498 temp = delsq2 / ( d( i )+sq2 )
500 work( j ) = d( j ) + d( i ) + temp
501 delta( j ) = ( d( j )-d( i ) ) - temp
506 psi = psi + z( j )*z( j ) / ( work( j )*delta( j ) )
510 DO 120 j = n, i + 2, -1
511 phi = phi + z( j )*z( j ) / ( work( j )*delta( j ) )
513 c = rhoinv + psi + phi
514 w = c + z( i )*z( i ) / ( work( i )*delta( i ) ) +
515 $ z( ip1 )*z( ip1 ) / ( work( ip1 )*delta( ip1 ) )
527 sgub = delsq2 / ( d( i )+sq2 )
528 a = c*delsq + z( i )*z( i ) + z( ip1 )*z( ip1 )
529 b = z( i )*z( i )*delsq
531 tau2 = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
533 tau2 = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
540 tau = tau2 / ( d( i )+sqrt( d( i )*d( i )+tau2 ) )
542 IF( (d(i).LE.temp*d(ip1)).AND.(abs(z(i)).LE.temp)
543 $ .AND.(d(i).GT.zero) )
THEN
544 tau = min( ten*d(i), sgub )
555 sglb = -delsq2 / ( d( ii )+sq2 )
557 a = c*delsq - z( i )*z( i ) - z( ip1 )*z( ip1 )
558 b = z( ip1 )*z( ip1 )*delsq
560 tau2 = two*b / ( a-sqrt( abs( a*a+four*b*c ) ) )
562 tau2 = -( a+sqrt( abs( a*a+four*b*c ) ) ) / ( two*c )
569 tau = tau2 / ( d( ip1 )+sqrt( abs( d( ip1 )*d( ip1 )+
573 sigma = d( ii ) + tau
575 work( j ) = d( j ) + d( ii ) + tau
576 delta( j ) = ( d( j )-d( ii ) ) - tau
587 temp = z( j ) / ( work( j )*delta( j ) )
588 psi = psi + z( j )*temp
589 dpsi = dpsi + temp*temp
590 erretm = erretm + psi
592 erretm = abs( erretm )
598 DO 160 j = n, iip1, -1
599 temp = z( j ) / ( work( j )*delta( j ) )
600 phi = phi + z( j )*temp
601 dphi = dphi + temp*temp
602 erretm = erretm + phi
605 w = rhoinv + phi + psi
618 IF( ii.EQ.1 .OR. ii.EQ.n )
621 temp = z( ii ) / ( work( ii )*delta( ii ) )
622 dw = dpsi + dphi + temp*temp
625 erretm = eight*( phi-psi ) + erretm + two*rhoinv
626 $ + three*abs( temp )
631 IF( abs( w ).LE.eps*erretm )
THEN
636 sglb = max( sglb, tau )
638 sgub = min( sgub, tau )
644 IF( .NOT.swtch3 )
THEN
645 dtipsq = work( ip1 )*delta( ip1 )
646 dtisq = work( i )*delta( i )
648 c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
650 c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
652 a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
657 a = z( i )*z( i ) + dtipsq*dtipsq*( dpsi+dphi )
659 a = z( ip1 )*z( ip1 ) + dtisq*dtisq*( dpsi+dphi )
663 ELSE IF( a.LE.zero )
THEN
664 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
666 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
672 dtiim = work( iim1 )*delta( iim1 )
673 dtiip = work( iip1 )*delta( iip1 )
674 temp = rhoinv + psi + phi
676 temp1 = z( iim1 ) / dtiim
678 c = ( temp - dtiip*( dpsi+dphi ) ) -
679 $ ( d( iim1 )-d( iip1 ) )*( d( iim1 )+d( iip1 ) )*temp1
680 zz( 1 ) = z( iim1 )*z( iim1 )
681 IF( dpsi.LT.temp1 )
THEN
682 zz( 3 ) = dtiip*dtiip*dphi
684 zz( 3 ) = dtiip*dtiip*( ( dpsi-temp1 )+dphi )
687 temp1 = z( iip1 ) / dtiip
689 c = ( temp - dtiim*( dpsi+dphi ) ) -
690 $ ( d( iip1 )-d( iim1 ) )*( d( iim1 )+d( iip1 ) )*temp1
691 IF( dphi.LT.temp1 )
THEN
692 zz( 1 ) = dtiim*dtiim*dpsi
694 zz( 1 ) = dtiim*dtiim*( dpsi+( dphi-temp1 ) )
696 zz( 3 ) = z( iip1 )*z( iip1 )
698 zz( 2 ) = z( ii )*z( ii )
700 dd( 2 ) = delta( ii )*work( ii )
702 CALL slaed6( niter, orgati, c, dd, zz, w, eta, info )
711 dtipsq = work( ip1 )*delta( ip1 )
712 dtisq = work( i )*delta( i )
714 c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
716 c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
718 a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
723 a = z( i )*z( i ) + dtipsq*dtipsq*( dpsi+dphi )
725 a = z( ip1 )*z( ip1 ) + dtisq*dtisq*( dpsi+dphi)
729 ELSE IF( a.LE.zero )
THEN
730 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
732 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
746 eta = eta / ( sigma+sqrt( sigma*sigma+eta ) )
748 IF( temp.GT.sgub .OR. temp.LT.sglb )
THEN
750 eta = ( sgub-tau ) / two
752 eta = ( sglb-tau ) / two
755 IF( w .LT. zero )
THEN
756 IF( tau .GT. zero )
THEN
757 eta = sqrt(sgub*tau)-tau
760 IF( sglb .GT. zero )
THEN
761 eta = sqrt(sglb*tau)-tau
773 work( j ) = work( j ) + eta
774 delta( j ) = delta( j ) - eta
783 temp = z( j ) / ( work( j )*delta( j ) )
784 psi = psi + z( j )*temp
785 dpsi = dpsi + temp*temp
786 erretm = erretm + psi
788 erretm = abs( erretm )
794 DO 190 j = n, iip1, -1
795 temp = z( j ) / ( work( j )*delta( j ) )
796 phi = phi + z( j )*temp
797 dphi = dphi + temp*temp
798 erretm = erretm + phi
801 tau2 = work( ii )*delta( ii )
802 temp = z( ii ) / tau2
803 dw = dpsi + dphi + temp*temp
805 w = rhoinv + phi + psi + temp
806 erretm = eight*( phi-psi ) + erretm + two*rhoinv
807 $ + three*abs( temp )
812 IF( -w.GT.abs( prew ) / ten )
815 IF( w.GT.abs( prew ) / ten )
823 DO 230 niter = iter, maxit
827 IF( abs( w ).LE.eps*erretm )
THEN
833 sglb = max( sglb, tau )
835 sgub = min( sgub, tau )
840 IF( .NOT.swtch3 )
THEN
841 dtipsq = work( ip1 )*delta( ip1 )
842 dtisq = work( i )*delta( i )
843 IF( .NOT.swtch )
THEN
845 c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
847 c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
850 temp = z( ii ) / ( work( ii )*delta( ii ) )
852 dpsi = dpsi + temp*temp
854 dphi = dphi + temp*temp
856 c = w - dtisq*dpsi - dtipsq*dphi
858 a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
862 IF( .NOT.swtch )
THEN
864 a = z( i )*z( i ) + dtipsq*dtipsq*
867 a = z( ip1 )*z( ip1 ) +
868 $ dtisq*dtisq*( dpsi+dphi )
871 a = dtisq*dtisq*dpsi + dtipsq*dtipsq*dphi
875 ELSE IF( a.LE.zero )
THEN
876 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
878 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
884 dtiim = work( iim1 )*delta( iim1 )
885 dtiip = work( iip1 )*delta( iip1 )
886 temp = rhoinv + psi + phi
888 c = temp - dtiim*dpsi - dtiip*dphi
889 zz( 1 ) = dtiim*dtiim*dpsi
890 zz( 3 ) = dtiip*dtiip*dphi
893 temp1 = z( iim1 ) / dtiim
895 temp2 = ( d( iim1 )-d( iip1 ) )*
896 $ ( d( iim1 )+d( iip1 ) )*temp1
897 c = temp - dtiip*( dpsi+dphi ) - temp2
898 zz( 1 ) = z( iim1 )*z( iim1 )
899 IF( dpsi.LT.temp1 )
THEN
900 zz( 3 ) = dtiip*dtiip*dphi
902 zz( 3 ) = dtiip*dtiip*( ( dpsi-temp1 )+dphi )
905 temp1 = z( iip1 ) / dtiip
907 temp2 = ( d( iip1 )-d( iim1 ) )*
908 $ ( d( iim1 )+d( iip1 ) )*temp1
909 c = temp - dtiim*( dpsi+dphi ) - temp2
910 IF( dphi.LT.temp1 )
THEN
911 zz( 1 ) = dtiim*dtiim*dpsi
913 zz( 1 ) = dtiim*dtiim*( dpsi+( dphi-temp1 ) )
915 zz( 3 ) = z( iip1 )*z( iip1 )
919 dd( 2 ) = delta( ii )*work( ii )
921 CALL slaed6( niter, orgati, c, dd, zz, w, eta, info )
930 dtipsq = work( ip1 )*delta( ip1 )
931 dtisq = work( i )*delta( i )
932 IF( .NOT.swtch )
THEN
934 c = w - dtipsq*dw + delsq*( z( i )/dtisq )**2
936 c = w - dtisq*dw - delsq*( z( ip1 )/dtipsq )**2
939 temp = z( ii ) / ( work( ii )*delta( ii ) )
941 dpsi = dpsi + temp*temp
943 dphi = dphi + temp*temp
945 c = w - dtisq*dpsi - dtipsq*dphi
947 a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
951 IF( .NOT.swtch )
THEN
953 a = z( i )*z( i ) + dtipsq*dtipsq*
956 a = z( ip1 )*z( ip1 ) +
957 $ dtisq*dtisq*( dpsi+dphi )
960 a = dtisq*dtisq*dpsi + dtipsq*dtipsq*dphi
964 ELSE IF( a.LE.zero )
THEN
965 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
967 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
981 eta = eta / ( sigma+sqrt( sigma*sigma+eta ) )
983 IF( temp.GT.sgub .OR. temp.LT.sglb )
THEN
985 eta = ( sgub-tau ) / two
987 eta = ( sglb-tau ) / two
990 IF( w .LT. zero )
THEN
991 IF( tau .GT. zero )
THEN
992 eta = sqrt(sgub*tau)-tau
995 IF( sglb .GT. zero )
THEN
996 eta = sqrt(sglb*tau)-tau
1008 work( j ) = work( j ) + eta
1009 delta( j ) = delta( j ) - eta
1018 temp = z( j ) / ( work( j )*delta( j ) )
1019 psi = psi + z( j )*temp
1020 dpsi = dpsi + temp*temp
1021 erretm = erretm + psi
1023 erretm = abs( erretm )
1029 DO 220 j = n, iip1, -1
1030 temp = z( j ) / ( work( j )*delta( j ) )
1031 phi = phi + z( j )*temp
1032 dphi = dphi + temp*temp
1033 erretm = erretm + phi
1036 tau2 = work( ii )*delta( ii )
1037 temp = z( ii ) / tau2
1038 dw = dpsi + dphi + temp*temp
1040 w = rhoinv + phi + psi + temp
1041 erretm = eight*( phi-psi ) + erretm + two*rhoinv
1042 $ + three*abs( temp )
1045 IF( w*prew.GT.zero .AND. abs( w ).GT.abs( prew ) / ten )
1046 $ swtch = .NOT.swtch
subroutine slasd5(I, D, Z, DELTA, RHO, DSIGMA, WORK)
SLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification ...
subroutine slaed6(KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)
SLAED6 used by sstedc. Computes one Newton step in solution of the secular equation.
subroutine slasd4(N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO)
SLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modif...