144 SUBROUTINE slaed4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
155 REAL D( * ), DELTA( * ), Z( * )
162 parameter( maxit = 30 )
163 REAL ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
164 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
165 $ three = 3.0e0, four = 4.0e0, eight = 8.0e0,
169 LOGICAL ORGATI, SWTCH, SWTCH3
170 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
171 REAL A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
172 $ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
173 $ RHOINV, TAU, TEMP, TEMP1, W
186 INTRINSIC abs, max, min, sqrt
200 dlam = d( 1 ) + rho*z( 1 )*z( 1 )
205 CALL slaed5( i, d, z, delta, rho, dlam )
211 eps = slamch(
'Epsilon' )
231 delta( j ) = ( d( j )-d( i ) ) - midpt
236 psi = psi + z( j )*z( j ) / delta( j )
240 w = c + z( ii )*z( ii ) / delta( ii ) +
241 $ z( n )*z( n ) / delta( n )
244 temp = z( n-1 )*z( n-1 ) / ( d( n )-d( n-1 )+rho ) +
245 $ z( n )*z( n ) / rho
249 del = d( n ) - d( n-1 )
250 a = -c*del + z( n-1 )*z( n-1 ) + z( n )*z( n )
251 b = z( n )*z( n )*del
253 tau = two*b / ( sqrt( a*a+four*b*c )-a )
255 tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
265 del = d( n ) - d( n-1 )
266 a = -c*del + z( n-1 )*z( n-1 ) + z( n )*z( n )
267 b = z( n )*z( n )*del
269 tau = two*b / ( sqrt( a*a+four*b*c )-a )
271 tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
282 delta( j ) = ( d( j )-d( i ) ) - tau
291 temp = z( j ) / delta( j )
292 psi = psi + z( j )*temp
293 dpsi = dpsi + temp*temp
294 erretm = erretm + psi
296 erretm = abs( erretm )
300 temp = z( n ) / delta( n )
303 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
304 $ abs( tau )*( dpsi+dphi )
306 w = rhoinv + phi + psi
310 IF( abs( w ).LE.eps*erretm )
THEN
316 dltlb = max( dltlb, tau )
318 dltub = min( dltub, tau )
324 c = w - delta( n-1 )*dpsi - delta( n )*dphi
325 a = ( delta( n-1 )+delta( n ) )*w -
326 $ delta( n-1 )*delta( n )*( dpsi+dphi )
327 b = delta( n-1 )*delta( n )*w
334 ELSE IF( a.GE.zero )
THEN
335 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
337 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
347 $ eta = -w / ( dpsi+dphi )
349 IF( temp.GT.dltub .OR. temp.LT.dltlb )
THEN
351 eta = ( dltub-tau ) / two
353 eta = ( dltlb-tau ) / two
357 delta( j ) = delta( j ) - eta
368 temp = z( j ) / delta( j )
369 psi = psi + z( j )*temp
370 dpsi = dpsi + temp*temp
371 erretm = erretm + psi
373 erretm = abs( erretm )
377 temp = z( n ) / delta( n )
380 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
381 $ abs( tau )*( dpsi+dphi )
383 w = rhoinv + phi + psi
389 DO 90 niter = iter, maxit
393 IF( abs( w ).LE.eps*erretm )
THEN
399 dltlb = max( dltlb, tau )
401 dltub = min( dltub, tau )
406 c = w - delta( n-1 )*dpsi - delta( n )*dphi
407 a = ( delta( n-1 )+delta( n ) )*w -
408 $ delta( n-1 )*delta( n )*( dpsi+dphi )
409 b = delta( n-1 )*delta( n )*w
411 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
413 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
423 $ eta = -w / ( dpsi+dphi )
425 IF( temp.GT.dltub .OR. temp.LT.dltlb )
THEN
427 eta = ( dltub-tau ) / two
429 eta = ( dltlb-tau ) / two
433 delta( j ) = delta( j ) - eta
444 temp = z( j ) / delta( j )
445 psi = psi + z( j )*temp
446 dpsi = dpsi + temp*temp
447 erretm = erretm + psi
449 erretm = abs( erretm )
453 temp = z( n ) / delta( n )
456 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
457 $ abs( tau )*( dpsi+dphi )
459 w = rhoinv + phi + psi
479 del = d( ip1 ) - d( i )
482 delta( j ) = ( d( j )-d( i ) ) - midpt
487 psi = psi + z( j )*z( j ) / delta( j )
491 DO 120 j = n, i + 2, -1
492 phi = phi + z( j )*z( j ) / delta( j )
494 c = rhoinv + psi + phi
495 w = c + z( i )*z( i ) / delta( i ) +
496 $ z( ip1 )*z( ip1 ) / delta( ip1 )
505 a = c*del + z( i )*z( i ) + z( ip1 )*z( ip1 )
506 b = z( i )*z( i )*del
508 tau = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
510 tau = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
521 a = c*del - z( i )*z( i ) - z( ip1 )*z( ip1 )
522 b = z( ip1 )*z( ip1 )*del
524 tau = two*b / ( a-sqrt( abs( a*a+four*b*c ) ) )
526 tau = -( a+sqrt( abs( a*a+four*b*c ) ) ) / ( two*c )
534 delta( j ) = ( d( j )-d( i ) ) - tau
538 delta( j ) = ( d( j )-d( ip1 ) ) - tau
555 temp = z( j ) / delta( j )
556 psi = psi + z( j )*temp
557 dpsi = dpsi + temp*temp
558 erretm = erretm + psi
560 erretm = abs( erretm )
566 DO 160 j = n, iip1, -1
567 temp = z( j ) / delta( j )
568 phi = phi + z( j )*temp
569 dphi = dphi + temp*temp
570 erretm = erretm + phi
573 w = rhoinv + phi + psi
586 IF( ii.EQ.1 .OR. ii.EQ.n )
589 temp = z( ii ) / delta( ii )
590 dw = dpsi + dphi + temp*temp
593 erretm = eight*( phi-psi ) + erretm + two*rhoinv +
594 $ three*abs( temp ) + abs( tau )*dw
598 IF( abs( w ).LE.eps*erretm )
THEN
602 dlam = d( ip1 ) + tau
608 dltlb = max( dltlb, tau )
610 dltub = min( dltub, tau )
616 IF( .NOT.swtch3 )
THEN
618 c = w - delta( ip1 )*dw - ( d( i )-d( ip1 ) )*
619 $ ( z( i ) / delta( i ) )**2
621 c = w - delta( i )*dw - ( d( ip1 )-d( i ) )*
622 $ ( z( ip1 ) / delta( ip1 ) )**2
624 a = ( delta( i )+delta( ip1 ) )*w -
625 $ delta( i )*delta( ip1 )*dw
626 b = delta( i )*delta( ip1 )*w
630 a = z( i )*z( i ) + delta( ip1 )*delta( ip1 )*
633 a = z( ip1 )*z( ip1 ) + delta( i )*delta( i )*
638 ELSE IF( a.LE.zero )
THEN
639 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
641 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
647 temp = rhoinv + psi + phi
649 temp1 = z( iim1 ) / delta( iim1 )
651 c = temp - delta( iip1 )*( dpsi+dphi ) -
652 $ ( d( iim1 )-d( iip1 ) )*temp1
653 zz( 1 ) = z( iim1 )*z( iim1 )
654 zz( 3 ) = delta( iip1 )*delta( iip1 )*
655 $ ( ( dpsi-temp1 )+dphi )
657 temp1 = z( iip1 ) / delta( iip1 )
659 c = temp - delta( iim1 )*( dpsi+dphi ) -
660 $ ( d( iip1 )-d( iim1 ) )*temp1
661 zz( 1 ) = delta( iim1 )*delta( iim1 )*
662 $ ( dpsi+( dphi-temp1 ) )
663 zz( 3 ) = z( iip1 )*z( iip1 )
665 zz( 2 ) = z( ii )*z( ii )
666 CALL slaed6( niter, orgati, c, delta( iim1 ), zz, w, eta,
681 IF( temp.GT.dltub .OR. temp.LT.dltlb )
THEN
683 eta = ( dltub-tau ) / two
685 eta = ( dltlb-tau ) / two
692 delta( j ) = delta( j ) - eta
701 temp = z( j ) / delta( j )
702 psi = psi + z( j )*temp
703 dpsi = dpsi + temp*temp
704 erretm = erretm + psi
706 erretm = abs( erretm )
712 DO 200 j = n, iip1, -1
713 temp = z( j ) / delta( j )
714 phi = phi + z( j )*temp
715 dphi = dphi + temp*temp
716 erretm = erretm + phi
719 temp = z( ii ) / delta( ii )
720 dw = dpsi + dphi + temp*temp
722 w = rhoinv + phi + psi + temp
723 erretm = eight*( phi-psi ) + erretm + two*rhoinv +
724 $ three*abs( temp ) + abs( tau+eta )*dw
728 IF( -w.GT.abs( prew ) / ten )
731 IF( w.GT.abs( prew ) / ten )
741 DO 240 niter = iter, maxit
745 IF( abs( w ).LE.eps*erretm )
THEN
749 dlam = d( ip1 ) + tau
755 dltlb = max( dltlb, tau )
757 dltub = min( dltub, tau )
762 IF( .NOT.swtch3 )
THEN
763 IF( .NOT.swtch )
THEN
765 c = w - delta( ip1 )*dw -
766 $ ( d( i )-d( ip1 ) )*( z( i ) / delta( i ) )**2
768 c = w - delta( i )*dw - ( d( ip1 )-d( i ) )*
769 $ ( z( ip1 ) / delta( ip1 ) )**2
772 temp = z( ii ) / delta( ii )
774 dpsi = dpsi + temp*temp
776 dphi = dphi + temp*temp
778 c = w - delta( i )*dpsi - delta( ip1 )*dphi
780 a = ( delta( i )+delta( ip1 ) )*w -
781 $ delta( i )*delta( ip1 )*dw
782 b = delta( i )*delta( ip1 )*w
785 IF( .NOT.swtch )
THEN
787 a = z( i )*z( i ) + delta( ip1 )*
788 $ delta( ip1 )*( dpsi+dphi )
790 a = z( ip1 )*z( ip1 ) +
791 $ delta( i )*delta( i )*( dpsi+dphi )
794 a = delta( i )*delta( i )*dpsi +
795 $ delta( ip1 )*delta( ip1 )*dphi
799 ELSE IF( a.LE.zero )
THEN
800 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
802 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
808 temp = rhoinv + psi + phi
810 c = temp - delta( iim1 )*dpsi - delta( iip1 )*dphi
811 zz( 1 ) = delta( iim1 )*delta( iim1 )*dpsi
812 zz( 3 ) = delta( iip1 )*delta( iip1 )*dphi
815 temp1 = z( iim1 ) / delta( iim1 )
817 c = temp - delta( iip1 )*( dpsi+dphi ) -
818 $ ( d( iim1 )-d( iip1 ) )*temp1
819 zz( 1 ) = z( iim1 )*z( iim1 )
820 zz( 3 ) = delta( iip1 )*delta( iip1 )*
821 $ ( ( dpsi-temp1 )+dphi )
823 temp1 = z( iip1 ) / delta( iip1 )
825 c = temp - delta( iim1 )*( dpsi+dphi ) -
826 $ ( d( iip1 )-d( iim1 ) )*temp1
827 zz( 1 ) = delta( iim1 )*delta( iim1 )*
828 $ ( dpsi+( dphi-temp1 ) )
829 zz( 3 ) = z( iip1 )*z( iip1 )
832 CALL slaed6( niter, orgati, c, delta( iim1 ), zz, w, eta,
847 IF( temp.GT.dltub .OR. temp.LT.dltlb )
THEN
849 eta = ( dltub-tau ) / two
851 eta = ( dltlb-tau ) / two
856 delta( j ) = delta( j ) - eta
868 temp = z( j ) / delta( j )
869 psi = psi + z( j )*temp
870 dpsi = dpsi + temp*temp
871 erretm = erretm + psi
873 erretm = abs( erretm )
879 DO 230 j = n, iip1, -1
880 temp = z( j ) / delta( j )
881 phi = phi + z( j )*temp
882 dphi = dphi + temp*temp
883 erretm = erretm + phi
886 temp = z( ii ) / delta( ii )
887 dw = dpsi + dphi + temp*temp
889 w = rhoinv + phi + psi + temp
890 erretm = eight*( phi-psi ) + erretm + two*rhoinv +
891 $ three*abs( temp ) + abs( tau )*dw
892 IF( w*prew.GT.zero .AND. abs( w ).GT.abs( prew ) / ten )
903 dlam = d( ip1 ) + tau
subroutine slaed4(N, I, D, Z, DELTA, RHO, DLAM, INFO)
SLAED4 used by SSTEDC. Finds a single root of the secular equation.
subroutine slaed5(I, D, Z, DELTA, RHO, DLAM)
SLAED5 used by SSTEDC. Solves the 2-by-2 secular equation.
subroutine slaed6(KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)
SLAED6 used by SSTEDC. Computes one Newton step in solution of the secular equation.