158 REAL d( * ), delta( * ), z( * )
165 parameter ( maxit = 30 )
166 REAL zero, one, two, three, four, eight, ten
167 parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
168 $ three = 3.0e0, four = 4.0e0, eight = 8.0e0,
172 LOGICAL orgati, swtch, swtch3
173 INTEGER ii, iim1, iip1, ip1, iter, j, niter
174 REAL a, b, c, del, dltlb, dltub, dphi, dpsi, dw,
175 $ eps, erretm, eta, midpt, phi, prew, psi,
176 $ rhoinv, tau, temp, temp1, w
189 INTRINSIC abs, max, min, sqrt
203 dlam = d( 1 ) + rho*z( 1 )*z( 1 )
208 CALL slaed5( i, d, z, delta, rho, dlam )
234 delta( j ) = ( d( j )-d( i ) ) - midpt
239 psi = psi + z( j )*z( j ) / delta( j )
243 w = c + z( ii )*z( ii ) / delta( ii ) +
244 $ z( n )*z( n ) / delta( n )
247 temp = z( n-1 )*z( n-1 ) / ( d( n )-d( n-1 )+rho ) +
248 $ z( n )*z( n ) / rho
252 del = d( n ) - d( n-1 )
253 a = -c*del + z( n-1 )*z( n-1 ) + z( n )*z( n )
254 b = z( n )*z( n )*del
256 tau = two*b / ( sqrt( a*a+four*b*c )-a )
258 tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
268 del = d( n ) - d( n-1 )
269 a = -c*del + z( n-1 )*z( n-1 ) + z( n )*z( n )
270 b = z( n )*z( n )*del
272 tau = two*b / ( sqrt( a*a+four*b*c )-a )
274 tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
285 delta( j ) = ( d( j )-d( i ) ) - tau
294 temp = z( j ) / delta( j )
295 psi = psi + z( j )*temp
296 dpsi = dpsi + temp*temp
297 erretm = erretm + psi
299 erretm = abs( erretm )
303 temp = z( n ) / delta( n )
306 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
307 $ abs( tau )*( dpsi+dphi )
309 w = rhoinv + phi + psi
313 IF( abs( w ).LE.eps*erretm )
THEN
319 dltlb = max( dltlb, tau )
321 dltub = min( dltub, tau )
327 c = w - delta( n-1 )*dpsi - delta( n )*dphi
328 a = ( delta( n-1 )+delta( n ) )*w -
329 $ delta( n-1 )*delta( n )*( dpsi+dphi )
330 b = delta( n-1 )*delta( n )*w
337 ELSE IF( a.GE.zero )
THEN
338 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
340 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
350 $ eta = -w / ( dpsi+dphi )
352 IF( temp.GT.dltub .OR. temp.LT.dltlb )
THEN
354 eta = ( dltub-tau ) / two
356 eta = ( dltlb-tau ) / two
360 delta( j ) = delta( j ) - eta
371 temp = z( j ) / delta( j )
372 psi = psi + z( j )*temp
373 dpsi = dpsi + temp*temp
374 erretm = erretm + psi
376 erretm = abs( erretm )
380 temp = z( n ) / delta( n )
383 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
384 $ abs( tau )*( dpsi+dphi )
386 w = rhoinv + phi + psi
392 DO 90 niter = iter, maxit
396 IF( abs( w ).LE.eps*erretm )
THEN
402 dltlb = max( dltlb, tau )
404 dltub = min( dltub, tau )
409 c = w - delta( n-1 )*dpsi - delta( n )*dphi
410 a = ( delta( n-1 )+delta( n ) )*w -
411 $ delta( n-1 )*delta( n )*( dpsi+dphi )
412 b = delta( n-1 )*delta( n )*w
414 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
416 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
426 $ eta = -w / ( dpsi+dphi )
428 IF( temp.GT.dltub .OR. temp.LT.dltlb )
THEN
430 eta = ( dltub-tau ) / two
432 eta = ( dltlb-tau ) / two
436 delta( j ) = delta( j ) - eta
447 temp = z( j ) / delta( j )
448 psi = psi + z( j )*temp
449 dpsi = dpsi + temp*temp
450 erretm = erretm + psi
452 erretm = abs( erretm )
456 temp = z( n ) / delta( n )
459 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
460 $ abs( tau )*( dpsi+dphi )
462 w = rhoinv + phi + psi
482 del = d( ip1 ) - d( i )
485 delta( j ) = ( d( j )-d( i ) ) - midpt
490 psi = psi + z( j )*z( j ) / delta( j )
494 DO 120 j = n, i + 2, -1
495 phi = phi + z( j )*z( j ) / delta( j )
497 c = rhoinv + psi + phi
498 w = c + z( i )*z( i ) / delta( i ) +
499 $ z( ip1 )*z( ip1 ) / delta( ip1 )
508 a = c*del + z( i )*z( i ) + z( ip1 )*z( ip1 )
509 b = z( i )*z( i )*del
511 tau = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
513 tau = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
524 a = c*del - z( i )*z( i ) - z( ip1 )*z( ip1 )
525 b = z( ip1 )*z( ip1 )*del
527 tau = two*b / ( a-sqrt( abs( a*a+four*b*c ) ) )
529 tau = -( a+sqrt( abs( a*a+four*b*c ) ) ) / ( two*c )
537 delta( j ) = ( d( j )-d( i ) ) - tau
541 delta( j ) = ( d( j )-d( ip1 ) ) - tau
558 temp = z( j ) / delta( j )
559 psi = psi + z( j )*temp
560 dpsi = dpsi + temp*temp
561 erretm = erretm + psi
563 erretm = abs( erretm )
569 DO 160 j = n, iip1, -1
570 temp = z( j ) / delta( j )
571 phi = phi + z( j )*temp
572 dphi = dphi + temp*temp
573 erretm = erretm + phi
576 w = rhoinv + phi + psi
589 IF( ii.EQ.1 .OR. ii.EQ.n )
592 temp = z( ii ) / delta( ii )
593 dw = dpsi + dphi + temp*temp
596 erretm = eight*( phi-psi ) + erretm + two*rhoinv +
597 $ three*abs( temp ) + abs( tau )*dw
601 IF( abs( w ).LE.eps*erretm )
THEN
605 dlam = d( ip1 ) + tau
611 dltlb = max( dltlb, tau )
613 dltub = min( dltub, tau )
619 IF( .NOT.swtch3 )
THEN
621 c = w - delta( ip1 )*dw - ( d( i )-d( ip1 ) )*
622 $ ( z( i ) / delta( i ) )**2
624 c = w - delta( i )*dw - ( d( ip1 )-d( i ) )*
625 $ ( z( ip1 ) / delta( ip1 ) )**2
627 a = ( delta( i )+delta( ip1 ) )*w -
628 $ delta( i )*delta( ip1 )*dw
629 b = delta( i )*delta( ip1 )*w
633 a = z( i )*z( i ) + delta( ip1 )*delta( ip1 )*
636 a = z( ip1 )*z( ip1 ) + delta( i )*delta( i )*
641 ELSE IF( a.LE.zero )
THEN
642 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
644 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
650 temp = rhoinv + psi + phi
652 temp1 = z( iim1 ) / delta( iim1 )
654 c = temp - delta( iip1 )*( dpsi+dphi ) -
655 $ ( d( iim1 )-d( iip1 ) )*temp1
656 zz( 1 ) = z( iim1 )*z( iim1 )
657 zz( 3 ) = delta( iip1 )*delta( iip1 )*
658 $ ( ( dpsi-temp1 )+dphi )
660 temp1 = z( iip1 ) / delta( iip1 )
662 c = temp - delta( iim1 )*( dpsi+dphi ) -
663 $ ( d( iip1 )-d( iim1 ) )*temp1
664 zz( 1 ) = delta( iim1 )*delta( iim1 )*
665 $ ( dpsi+( dphi-temp1 ) )
666 zz( 3 ) = z( iip1 )*z( iip1 )
668 zz( 2 ) = z( ii )*z( ii )
669 CALL slaed6( niter, orgati, c, delta( iim1 ), zz, w, eta,
684 IF( temp.GT.dltub .OR. temp.LT.dltlb )
THEN
686 eta = ( dltub-tau ) / two
688 eta = ( dltlb-tau ) / two
695 delta( j ) = delta( j ) - eta
704 temp = z( j ) / delta( j )
705 psi = psi + z( j )*temp
706 dpsi = dpsi + temp*temp
707 erretm = erretm + psi
709 erretm = abs( erretm )
715 DO 200 j = n, iip1, -1
716 temp = z( j ) / delta( j )
717 phi = phi + z( j )*temp
718 dphi = dphi + temp*temp
719 erretm = erretm + phi
722 temp = z( ii ) / delta( ii )
723 dw = dpsi + dphi + temp*temp
725 w = rhoinv + phi + psi + temp
726 erretm = eight*( phi-psi ) + erretm + two*rhoinv +
727 $ three*abs( temp ) + abs( tau+eta )*dw
731 IF( -w.GT.abs( prew ) / ten )
734 IF( w.GT.abs( prew ) / ten )
744 DO 240 niter = iter, maxit
748 IF( abs( w ).LE.eps*erretm )
THEN
752 dlam = d( ip1 ) + tau
758 dltlb = max( dltlb, tau )
760 dltub = min( dltub, tau )
765 IF( .NOT.swtch3 )
THEN
766 IF( .NOT.swtch )
THEN
768 c = w - delta( ip1 )*dw -
769 $ ( d( i )-d( ip1 ) )*( z( i ) / delta( i ) )**2
771 c = w - delta( i )*dw - ( d( ip1 )-d( i ) )*
772 $ ( z( ip1 ) / delta( ip1 ) )**2
775 temp = z( ii ) / delta( ii )
777 dpsi = dpsi + temp*temp
779 dphi = dphi + temp*temp
781 c = w - delta( i )*dpsi - delta( ip1 )*dphi
783 a = ( delta( i )+delta( ip1 ) )*w -
784 $ delta( i )*delta( ip1 )*dw
785 b = delta( i )*delta( ip1 )*w
788 IF( .NOT.swtch )
THEN
790 a = z( i )*z( i ) + delta( ip1 )*
791 $ delta( ip1 )*( dpsi+dphi )
793 a = z( ip1 )*z( ip1 ) +
794 $ delta( i )*delta( i )*( dpsi+dphi )
797 a = delta( i )*delta( i )*dpsi +
798 $ delta( ip1 )*delta( ip1 )*dphi
802 ELSE IF( a.LE.zero )
THEN
803 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
805 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
811 temp = rhoinv + psi + phi
813 c = temp - delta( iim1 )*dpsi - delta( iip1 )*dphi
814 zz( 1 ) = delta( iim1 )*delta( iim1 )*dpsi
815 zz( 3 ) = delta( iip1 )*delta( iip1 )*dphi
818 temp1 = z( iim1 ) / delta( iim1 )
820 c = temp - delta( iip1 )*( dpsi+dphi ) -
821 $ ( d( iim1 )-d( iip1 ) )*temp1
822 zz( 1 ) = z( iim1 )*z( iim1 )
823 zz( 3 ) = delta( iip1 )*delta( iip1 )*
824 $ ( ( dpsi-temp1 )+dphi )
826 temp1 = z( iip1 ) / delta( iip1 )
828 c = temp - delta( iim1 )*( dpsi+dphi ) -
829 $ ( d( iip1 )-d( iim1 ) )*temp1
830 zz( 1 ) = delta( iim1 )*delta( iim1 )*
831 $ ( dpsi+( dphi-temp1 ) )
832 zz( 3 ) = z( iip1 )*z( iip1 )
835 CALL slaed6( niter, orgati, c, delta( iim1 ), zz, w, eta,
850 IF( temp.GT.dltub .OR. temp.LT.dltlb )
THEN
852 eta = ( dltub-tau ) / two
854 eta = ( dltlb-tau ) / two
859 delta( j ) = delta( j ) - eta
871 temp = z( j ) / delta( j )
872 psi = psi + z( j )*temp
873 dpsi = dpsi + temp*temp
874 erretm = erretm + psi
876 erretm = abs( erretm )
882 DO 230 j = n, iip1, -1
883 temp = z( j ) / delta( j )
884 phi = phi + z( j )*temp
885 dphi = dphi + temp*temp
886 erretm = erretm + phi
889 temp = z( ii ) / delta( ii )
890 dw = dpsi + dphi + temp*temp
892 w = rhoinv + phi + psi + temp
893 erretm = eight*( phi-psi ) + erretm + two*rhoinv +
894 $ three*abs( temp ) + abs( tau )*dw
895 IF( w*prew.GT.zero .AND. abs( w ).GT.abs( prew ) / ten )
906 dlam = d( ip1 ) + tau
subroutine slaed5(I, D, Z, DELTA, RHO, DLAM)
SLAED5 used by sstedc. Solves the 2-by-2 secular equation.
real function slamch(CMACH)
SLAMCH
subroutine slaed6(KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)
SLAED6 used by sstedc. Computes one Newton step in solution of the secular equation.