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 )
243 temp1 = temp / ( d( n )+sqrt( d( n )*d( n )+temp ) )
245 work( j ) = d( j ) + d( n ) + temp1
246 delta( j ) = ( d( j )-d( n ) ) - temp1
251 psi = psi + z( j )*z( j ) / ( delta( j )*work( j ) )
255 w = c + z( ii )*z( ii ) / ( delta( ii )*work( ii ) ) +
256 $ z( n )*z( n ) / ( delta( n )*work( n ) )
259 temp1 = sqrt( d( n )*d( n )+rho )
260 temp = z( n-1 )*z( n-1 ) / ( ( d( n-1 )+temp1 )*
261 $ ( d( n )-d( n-1 )+rho / ( d( n )+temp1 ) ) ) +
262 $ z( n )*z( n ) / rho
270 delsq = ( d( n )-d( n-1 ) )*( d( n )+d( n-1 ) )
271 a = -c*delsq + z( n-1 )*z( n-1 ) + z( n )*z( n )
272 b = z( n )*z( n )*delsq
274 tau2 = two*b / ( sqrt( a*a+four*b*c )-a )
276 tau2 = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
284 delsq = ( d( n )-d( n-1 ) )*( d( n )+d( n-1 ) )
285 a = -c*delsq + z( n-1 )*z( n-1 ) + z( n )*z( n )
286 b = z( n )*z( n )*delsq
292 tau2 = two*b / ( sqrt( a*a+four*b*c )-a )
294 tau2 = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
304 tau = tau2 / ( d( n )+sqrt( d( n )*d( n )+tau2 ) )
308 delta( j ) = ( d( j )-d( n ) ) - tau
309 work( j ) = d( j ) + d( n ) + tau
318 temp = z( j ) / ( delta( j )*work( j ) )
319 psi = psi + z( j )*temp
320 dpsi = dpsi + temp*temp
321 erretm = erretm + psi
323 erretm = abs( erretm )
327 temp = z( n ) / ( delta( n )*work( n ) )
330 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv
333 w = rhoinv + phi + psi
337 IF( abs( w ).LE.eps*erretm )
THEN
344 dtnsq1 = work( n-1 )*delta( n-1 )
345 dtnsq = work( n )*delta( n )
346 c = w - dtnsq1*dpsi - dtnsq*dphi
347 a = ( dtnsq+dtnsq1 )*w - dtnsq*dtnsq1*( dpsi+dphi )
352 eta = rho - sigma*sigma
353 ELSE IF( a.GE.zero )
THEN
354 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
356 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
366 $ eta = -w / ( dpsi+dphi )
371 eta = eta / ( sigma+sqrt( eta+sigma*sigma ) )
376 delta( j ) = delta( j ) - eta
377 work( j ) = work( j ) + eta
386 temp = z( j ) / ( work( j )*delta( j ) )
387 psi = psi + z( j )*temp
388 dpsi = dpsi + temp*temp
389 erretm = erretm + psi
391 erretm = abs( erretm )
395 tau2 = work( n )*delta( n )
399 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv
402 w = rhoinv + phi + psi
408 DO 90 niter = iter, maxit
412 IF( abs( w ).LE.eps*erretm )
THEN
418 dtnsq1 = work( n-1 )*delta( n-1 )
419 dtnsq = work( n )*delta( n )
420 c = w - dtnsq1*dpsi - dtnsq*dphi
421 a = ( dtnsq+dtnsq1 )*w - dtnsq1*dtnsq*( dpsi+dphi )
424 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
426 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
436 $ eta = -w / ( dpsi+dphi )
441 eta = eta / ( sigma+sqrt( eta+sigma*sigma ) )
446 delta( j ) = delta( j ) - eta
447 work( j ) = work( j ) + eta
456 temp = z( j ) / ( work( j )*delta( j ) )
457 psi = psi + z( j )*temp
458 dpsi = dpsi + temp*temp
459 erretm = erretm + psi
461 erretm = abs( erretm )
465 tau2 = work( n )*delta( n )
469 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv
472 w = rhoinv + phi + psi
491 delsq = ( d( ip1 )-d( i ) )*( d( ip1 )+d( i ) )
493 sq2=sqrt( ( d( i )*d( i )+d( ip1 )*d( ip1 ) ) / two )
494 temp = delsq2 / ( d( i )+sq2 )
496 work( j ) = d( j ) + d( i ) + temp
497 delta( j ) = ( d( j )-d( i ) ) - temp
502 psi = psi + z( j )*z( j ) / ( work( j )*delta( j ) )
506 DO 120 j = n, i + 2, -1
507 phi = phi + z( j )*z( j ) / ( work( j )*delta( j ) )
509 c = rhoinv + psi + phi
510 w = c + z( i )*z( i ) / ( work( i )*delta( i ) ) +
511 $ z( ip1 )*z( ip1 ) / ( work( ip1 )*delta( ip1 ) )
523 sgub = delsq2 / ( d( i )+sq2 )
524 a = c*delsq + z( i )*z( i ) + z( ip1 )*z( ip1 )
525 b = z( i )*z( i )*delsq
527 tau2 = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
529 tau2 = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
536 tau = tau2 / ( d( i )+sqrt( d( i )*d( i )+tau2 ) )
538 IF( (d(i).LE.temp*d(ip1)).AND.(abs(z(i)).LE.temp)
539 $ .AND.(d(i).GT.zero) )
THEN
540 tau = min( ten*d(i), sgub )
551 sglb = -delsq2 / ( d( ii )+sq2 )
553 a = c*delsq - z( i )*z( i ) - z( ip1 )*z( ip1 )
554 b = z( ip1 )*z( ip1 )*delsq
556 tau2 = two*b / ( a-sqrt( abs( a*a+four*b*c ) ) )
558 tau2 = -( a+sqrt( abs( a*a+four*b*c ) ) ) / ( two*c )
565 tau = tau2 / ( d( ip1 )+sqrt( abs( d( ip1 )*d( ip1 )+
569 sigma = d( ii ) + tau
571 work( j ) = d( j ) + d( ii ) + tau
572 delta( j ) = ( d( j )-d( ii ) ) - tau
583 temp = z( j ) / ( work( j )*delta( j ) )
584 psi = psi + z( j )*temp
585 dpsi = dpsi + temp*temp
586 erretm = erretm + psi
588 erretm = abs( erretm )
594 DO 160 j = n, iip1, -1
595 temp = z( j ) / ( work( j )*delta( j ) )
596 phi = phi + z( j )*temp
597 dphi = dphi + temp*temp
598 erretm = erretm + phi
601 w = rhoinv + phi + psi
614 IF( ii.EQ.1 .OR. ii.EQ.n )
617 temp = z( ii ) / ( work( ii )*delta( ii ) )
618 dw = dpsi + dphi + temp*temp
621 erretm = eight*( phi-psi ) + erretm + two*rhoinv
622 $ + three*abs( temp )
627 IF( abs( w ).LE.eps*erretm )
THEN
632 sglb = max( sglb, tau )
634 sgub = min( sgub, tau )
640 IF( .NOT.swtch3 )
THEN
641 dtipsq = work( ip1 )*delta( ip1 )
642 dtisq = work( i )*delta( i )
644 c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
646 c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
648 a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
653 a = z( i )*z( i ) + dtipsq*dtipsq*( dpsi+dphi )
655 a = z( ip1 )*z( ip1 ) + dtisq*dtisq*( dpsi+dphi )
659 ELSE IF( a.LE.zero )
THEN
660 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
662 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
668 dtiim = work( iim1 )*delta( iim1 )
669 dtiip = work( iip1 )*delta( iip1 )
670 temp = rhoinv + psi + phi
672 temp1 = z( iim1 ) / dtiim
674 c = ( temp - dtiip*( dpsi+dphi ) ) -
675 $ ( d( iim1 )-d( iip1 ) )*( d( iim1 )+d( iip1 ) )*temp1
676 zz( 1 ) = z( iim1 )*z( iim1 )
677 IF( dpsi.LT.temp1 )
THEN
678 zz( 3 ) = dtiip*dtiip*dphi
680 zz( 3 ) = dtiip*dtiip*( ( dpsi-temp1 )+dphi )
683 temp1 = z( iip1 ) / dtiip
685 c = ( temp - dtiim*( dpsi+dphi ) ) -
686 $ ( d( iip1 )-d( iim1 ) )*( d( iim1 )+d( iip1 ) )*temp1
687 IF( dphi.LT.temp1 )
THEN
688 zz( 1 ) = dtiim*dtiim*dpsi
690 zz( 1 ) = dtiim*dtiim*( dpsi+( dphi-temp1 ) )
692 zz( 3 ) = z( iip1 )*z( iip1 )
694 zz( 2 ) = z( ii )*z( ii )
696 dd( 2 ) = delta( ii )*work( ii )
698 CALL
slaed6( niter, orgati, c, dd, zz, w, eta, info )
707 dtipsq = work( ip1 )*delta( ip1 )
708 dtisq = work( i )*delta( i )
710 c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
712 c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
714 a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
719 a = z( i )*z( i ) + dtipsq*dtipsq*( dpsi+dphi )
721 a = z( ip1 )*z( ip1 ) + dtisq*dtisq*( dpsi+dphi)
725 ELSE IF( a.LE.zero )
THEN
726 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
728 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
742 eta = eta / ( sigma+sqrt( sigma*sigma+eta ) )
744 IF( temp.GT.sgub .OR. temp.LT.sglb )
THEN
746 eta = ( sgub-tau ) / two
748 eta = ( sglb-tau ) / two
751 IF( w .LT. zero )
THEN
752 IF( tau .GT. zero )
THEN
753 eta = sqrt(sgub*tau)-tau
756 IF( sglb .GT. zero )
THEN
757 eta = sqrt(sglb*tau)-tau
769 work( j ) = work( j ) + eta
770 delta( j ) = delta( j ) - eta
779 temp = z( j ) / ( work( j )*delta( j ) )
780 psi = psi + z( j )*temp
781 dpsi = dpsi + temp*temp
782 erretm = erretm + psi
784 erretm = abs( erretm )
790 DO 190 j = n, iip1, -1
791 temp = z( j ) / ( work( j )*delta( j ) )
792 phi = phi + z( j )*temp
793 dphi = dphi + temp*temp
794 erretm = erretm + phi
797 tau2 = work( ii )*delta( ii )
798 temp = z( ii ) / tau2
799 dw = dpsi + dphi + temp*temp
801 w = rhoinv + phi + psi + temp
802 erretm = eight*( phi-psi ) + erretm + two*rhoinv
803 $ + three*abs( temp )
808 IF( -w.GT.abs( prew ) / ten )
811 IF( w.GT.abs( prew ) / ten )
819 DO 230 niter = iter, maxit
823 IF( abs( w ).LE.eps*erretm )
THEN
829 sglb = max( sglb, tau )
831 sgub = min( sgub, tau )
836 IF( .NOT.swtch3 )
THEN
837 dtipsq = work( ip1 )*delta( ip1 )
838 dtisq = work( i )*delta( i )
839 IF( .NOT.swtch )
THEN
841 c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
843 c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
846 temp = z( ii ) / ( work( ii )*delta( ii ) )
848 dpsi = dpsi + temp*temp
850 dphi = dphi + temp*temp
852 c = w - dtisq*dpsi - dtipsq*dphi
854 a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
858 IF( .NOT.swtch )
THEN
860 a = z( i )*z( i ) + dtipsq*dtipsq*
863 a = z( ip1 )*z( ip1 ) +
864 $ dtisq*dtisq*( dpsi+dphi )
867 a = dtisq*dtisq*dpsi + dtipsq*dtipsq*dphi
871 ELSE IF( a.LE.zero )
THEN
872 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
874 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
880 dtiim = work( iim1 )*delta( iim1 )
881 dtiip = work( iip1 )*delta( iip1 )
882 temp = rhoinv + psi + phi
884 c = temp - dtiim*dpsi - dtiip*dphi
885 zz( 1 ) = dtiim*dtiim*dpsi
886 zz( 3 ) = dtiip*dtiip*dphi
889 temp1 = z( iim1 ) / dtiim
891 temp2 = ( d( iim1 )-d( iip1 ) )*
892 $ ( d( iim1 )+d( iip1 ) )*temp1
893 c = temp - dtiip*( dpsi+dphi ) - temp2
894 zz( 1 ) = z( iim1 )*z( iim1 )
895 IF( dpsi.LT.temp1 )
THEN
896 zz( 3 ) = dtiip*dtiip*dphi
898 zz( 3 ) = dtiip*dtiip*( ( dpsi-temp1 )+dphi )
901 temp1 = z( iip1 ) / dtiip
903 temp2 = ( d( iip1 )-d( iim1 ) )*
904 $ ( d( iim1 )+d( iip1 ) )*temp1
905 c = temp - dtiim*( dpsi+dphi ) - temp2
906 IF( dphi.LT.temp1 )
THEN
907 zz( 1 ) = dtiim*dtiim*dpsi
909 zz( 1 ) = dtiim*dtiim*( dpsi+( dphi-temp1 ) )
911 zz( 3 ) = z( iip1 )*z( iip1 )
915 dd( 2 ) = delta( ii )*work( ii )
917 CALL
slaed6( niter, orgati, c, dd, zz, w, eta, info )
926 dtipsq = work( ip1 )*delta( ip1 )
927 dtisq = work( i )*delta( i )
928 IF( .NOT.swtch )
THEN
930 c = w - dtipsq*dw + delsq*( z( i )/dtisq )**2
932 c = w - dtisq*dw - delsq*( z( ip1 )/dtipsq )**2
935 temp = z( ii ) / ( work( ii )*delta( ii ) )
937 dpsi = dpsi + temp*temp
939 dphi = dphi + temp*temp
941 c = w - dtisq*dpsi - dtipsq*dphi
943 a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
947 IF( .NOT.swtch )
THEN
949 a = z( i )*z( i ) + dtipsq*dtipsq*
952 a = z( ip1 )*z( ip1 ) +
953 $ dtisq*dtisq*( dpsi+dphi )
956 a = dtisq*dtisq*dpsi + dtipsq*dtipsq*dphi
960 ELSE IF( a.LE.zero )
THEN
961 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
963 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
977 eta = eta / ( sigma+sqrt( sigma*sigma+eta ) )
979 IF( temp.GT.sgub .OR. temp.LT.sglb )
THEN
981 eta = ( sgub-tau ) / two
983 eta = ( sglb-tau ) / two
986 IF( w .LT. zero )
THEN
987 IF( tau .GT. zero )
THEN
988 eta = sqrt(sgub*tau)-tau
991 IF( sglb .GT. zero )
THEN
992 eta = sqrt(sglb*tau)-tau
1004 work( j ) = work( j ) + eta
1005 delta( j ) = delta( j ) - eta
1014 temp = z( j ) / ( work( j )*delta( j ) )
1015 psi = psi + z( j )*temp
1016 dpsi = dpsi + temp*temp
1017 erretm = erretm + psi
1019 erretm = abs( erretm )
1025 DO 220 j = n, iip1, -1
1026 temp = z( j ) / ( work( j )*delta( j ) )
1027 phi = phi + z( j )*temp
1028 dphi = dphi + temp*temp
1029 erretm = erretm + phi
1032 tau2 = work( ii )*delta( ii )
1033 temp = z( ii ) / tau2
1034 dw = dpsi + dphi + temp*temp
1036 w = rhoinv + phi + psi + temp
1037 erretm = eight*( phi-psi ) + erretm + two*rhoinv
1038 $ + three*abs( temp )
1041 IF( w*prew.GT.zero .AND. abs( w ).GT.abs( prew ) / ten )
1042 $ swtch = .NOT.swtch