LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine slaed4 ( integer  N,
integer  I,
real, dimension( * )  D,
real, dimension( * )  Z,
real, dimension( * )  DELTA,
real  RHO,
real  DLAM,
integer  INFO 
)

SLAED4 used by sstedc. Finds a single root of the secular equation.

Download SLAED4 + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 This subroutine computes the I-th updated eigenvalue of a symmetric
 rank-one modification to a diagonal matrix whose elements are
 given in the array d, and that

            D(i) < D(j)  for  i < j

 and that RHO > 0.  This is arranged by the calling routine, and is
 no loss in generality.  The rank-one modified system is thus

            diag( D )  +  RHO * Z * Z_transpose.

 where we assume the Euclidean norm of Z is 1.

 The method consists of approximating the rational functions in the
 secular equation by simpler interpolating rational functions.
Parameters
[in]N
          N is INTEGER
         The length of all arrays.
[in]I
          I is INTEGER
         The index of the eigenvalue to be computed.  1 <= I <= N.
[in]D
          D is REAL array, dimension (N)
         The original eigenvalues.  It is assumed that they are in
         order, D(I) < D(J)  for I < J.
[in]Z
          Z is REAL array, dimension (N)
         The components of the updating vector.
[out]DELTA
          DELTA is REAL array, dimension (N)
         If N .GT. 2, DELTA contains (D(j) - lambda_I) in its  j-th
         component.  If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5
         for detail. The vector DELTA contains the information necessary
         to construct the eigenvectors by SLAED3 and SLAED9.
[in]RHO
          RHO is REAL
         The scalar in the symmetric updating formula.
[out]DLAM
          DLAM is REAL
         The computed lambda_I, the I-th updated eigenvalue.
[out]INFO
          INFO is INTEGER
         = 0:  successful exit
         > 0:  if INFO = 1, the updating process failed.
Internal Parameters: 
  Logical variable ORGATI (origin-at-i?) is used for distinguishing
  whether D(i) or D(i+1) is treated as the origin.

            ORGATI = .true.    origin at i
            ORGATI = .false.   origin at i+1

   Logical variable SWTCH3 (switch-for-3-poles?) is for noting
   if we are working with THREE poles!

   MAXIT is the maximum number of iterations allowed for each
   eigenvalue.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

Definition at line 147 of file slaed4.f.

147 *
148 * -- LAPACK computational routine (version 3.4.2) --
149 * -- LAPACK is a software package provided by Univ. of Tennessee, --
150 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151 * September 2012
152 *
153 * .. Scalar Arguments ..
154  INTEGER i, info, n
155  REAL dlam, rho
156 * ..
157 * .. Array Arguments ..
158  REAL d( * ), delta( * ), z( * )
159 * ..
160 *
161 * =====================================================================
162 *
163 * .. Parameters ..
164  INTEGER maxit
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,
169  $ ten = 10.0e0 )
170 * ..
171 * .. Local Scalars ..
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
177 * ..
178 * .. Local Arrays ..
179  REAL zz( 3 )
180 * ..
181 * .. External Functions ..
182  REAL slamch
183  EXTERNAL slamch
184 * ..
185 * .. External Subroutines ..
186  EXTERNAL slaed5, slaed6
187 * ..
188 * .. Intrinsic Functions ..
189  INTRINSIC abs, max, min, sqrt
190 * ..
191 * .. Executable Statements ..
192 *
193 * Since this routine is called in an inner loop, we do no argument
194 * checking.
195 *
196 * Quick return for N=1 and 2.
197 *
198  info = 0
199  IF( n.EQ.1 ) THEN
200 *
201 * Presumably, I=1 upon entry
202 *
203  dlam = d( 1 ) + rho*z( 1 )*z( 1 )
204  delta( 1 ) = one
205  RETURN
206  END IF
207  IF( n.EQ.2 ) THEN
208  CALL slaed5( i, d, z, delta, rho, dlam )
209  RETURN
210  END IF
211 *
212 * Compute machine epsilon
213 *
214  eps = slamch( 'Epsilon' )
215  rhoinv = one / rho
216 *
217 * The case I = N
218 *
219  IF( i.EQ.n ) THEN
220 *
221 * Initialize some basic variables
222 *
223  ii = n - 1
224  niter = 1
225 *
226 * Calculate initial guess
227 *
228  midpt = rho / two
229 *
230 * If ||Z||_2 is not one, then TEMP should be set to
231 * RHO * ||Z||_2^2 / TWO
232 *
233  DO 10 j = 1, n
234  delta( j ) = ( d( j )-d( i ) ) - midpt
235  10 CONTINUE
236 *
237  psi = zero
238  DO 20 j = 1, n - 2
239  psi = psi + z( j )*z( j ) / delta( j )
240  20 CONTINUE
241 *
242  c = rhoinv + psi
243  w = c + z( ii )*z( ii ) / delta( ii ) +
244  $ z( n )*z( n ) / delta( n )
245 *
246  IF( w.LE.zero ) THEN
247  temp = z( n-1 )*z( n-1 ) / ( d( n )-d( n-1 )+rho ) +
248  $ z( n )*z( n ) / rho
249  IF( c.LE.temp ) THEN
250  tau = rho
251  ELSE
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
255  IF( a.LT.zero ) THEN
256  tau = two*b / ( sqrt( a*a+four*b*c )-a )
257  ELSE
258  tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
259  END IF
260  END IF
261 *
262 * It can be proved that
263 * D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO
264 *
265  dltlb = midpt
266  dltub = rho
267  ELSE
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
271  IF( a.LT.zero ) THEN
272  tau = two*b / ( sqrt( a*a+four*b*c )-a )
273  ELSE
274  tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
275  END IF
276 *
277 * It can be proved that
278 * D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2
279 *
280  dltlb = zero
281  dltub = midpt
282  END IF
283 *
284  DO 30 j = 1, n
285  delta( j ) = ( d( j )-d( i ) ) - tau
286  30 CONTINUE
287 *
288 * Evaluate PSI and the derivative DPSI
289 *
290  dpsi = zero
291  psi = zero
292  erretm = zero
293  DO 40 j = 1, ii
294  temp = z( j ) / delta( j )
295  psi = psi + z( j )*temp
296  dpsi = dpsi + temp*temp
297  erretm = erretm + psi
298  40 CONTINUE
299  erretm = abs( erretm )
300 *
301 * Evaluate PHI and the derivative DPHI
302 *
303  temp = z( n ) / delta( n )
304  phi = z( n )*temp
305  dphi = temp*temp
306  erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
307  $ abs( tau )*( dpsi+dphi )
308 *
309  w = rhoinv + phi + psi
310 *
311 * Test for convergence
312 *
313  IF( abs( w ).LE.eps*erretm ) THEN
314  dlam = d( i ) + tau
315  GO TO 250
316  END IF
317 *
318  IF( w.LE.zero ) THEN
319  dltlb = max( dltlb, tau )
320  ELSE
321  dltub = min( dltub, tau )
322  END IF
323 *
324 * Calculate the new step
325 *
326  niter = niter + 1
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
331  IF( c.LT.zero )
332  $ c = abs( c )
333  IF( c.EQ.zero ) THEN
334 * ETA = B/A
335 * ETA = RHO - TAU
336  eta = dltub - tau
337  ELSE IF( a.GE.zero ) THEN
338  eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
339  ELSE
340  eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
341  END IF
342 *
343 * Note, eta should be positive if w is negative, and
344 * eta should be negative otherwise. However,
345 * if for some reason caused by roundoff, eta*w > 0,
346 * we simply use one Newton step instead. This way
347 * will guarantee eta*w < 0.
348 *
349  IF( w*eta.GT.zero )
350  $ eta = -w / ( dpsi+dphi )
351  temp = tau + eta
352  IF( temp.GT.dltub .OR. temp.LT.dltlb ) THEN
353  IF( w.LT.zero ) THEN
354  eta = ( dltub-tau ) / two
355  ELSE
356  eta = ( dltlb-tau ) / two
357  END IF
358  END IF
359  DO 50 j = 1, n
360  delta( j ) = delta( j ) - eta
361  50 CONTINUE
362 *
363  tau = tau + eta
364 *
365 * Evaluate PSI and the derivative DPSI
366 *
367  dpsi = zero
368  psi = zero
369  erretm = zero
370  DO 60 j = 1, ii
371  temp = z( j ) / delta( j )
372  psi = psi + z( j )*temp
373  dpsi = dpsi + temp*temp
374  erretm = erretm + psi
375  60 CONTINUE
376  erretm = abs( erretm )
377 *
378 * Evaluate PHI and the derivative DPHI
379 *
380  temp = z( n ) / delta( n )
381  phi = z( n )*temp
382  dphi = temp*temp
383  erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
384  $ abs( tau )*( dpsi+dphi )
385 *
386  w = rhoinv + phi + psi
387 *
388 * Main loop to update the values of the array DELTA
389 *
390  iter = niter + 1
391 *
392  DO 90 niter = iter, maxit
393 *
394 * Test for convergence
395 *
396  IF( abs( w ).LE.eps*erretm ) THEN
397  dlam = d( i ) + tau
398  GO TO 250
399  END IF
400 *
401  IF( w.LE.zero ) THEN
402  dltlb = max( dltlb, tau )
403  ELSE
404  dltub = min( dltub, tau )
405  END IF
406 *
407 * Calculate the new step
408 *
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
413  IF( a.GE.zero ) THEN
414  eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
415  ELSE
416  eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
417  END IF
418 *
419 * Note, eta should be positive if w is negative, and
420 * eta should be negative otherwise. However,
421 * if for some reason caused by roundoff, eta*w > 0,
422 * we simply use one Newton step instead. This way
423 * will guarantee eta*w < 0.
424 *
425  IF( w*eta.GT.zero )
426  $ eta = -w / ( dpsi+dphi )
427  temp = tau + eta
428  IF( temp.GT.dltub .OR. temp.LT.dltlb ) THEN
429  IF( w.LT.zero ) THEN
430  eta = ( dltub-tau ) / two
431  ELSE
432  eta = ( dltlb-tau ) / two
433  END IF
434  END IF
435  DO 70 j = 1, n
436  delta( j ) = delta( j ) - eta
437  70 CONTINUE
438 *
439  tau = tau + eta
440 *
441 * Evaluate PSI and the derivative DPSI
442 *
443  dpsi = zero
444  psi = zero
445  erretm = zero
446  DO 80 j = 1, ii
447  temp = z( j ) / delta( j )
448  psi = psi + z( j )*temp
449  dpsi = dpsi + temp*temp
450  erretm = erretm + psi
451  80 CONTINUE
452  erretm = abs( erretm )
453 *
454 * Evaluate PHI and the derivative DPHI
455 *
456  temp = z( n ) / delta( n )
457  phi = z( n )*temp
458  dphi = temp*temp
459  erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
460  $ abs( tau )*( dpsi+dphi )
461 *
462  w = rhoinv + phi + psi
463  90 CONTINUE
464 *
465 * Return with INFO = 1, NITER = MAXIT and not converged
466 *
467  info = 1
468  dlam = d( i ) + tau
469  GO TO 250
470 *
471 * End for the case I = N
472 *
473  ELSE
474 *
475 * The case for I < N
476 *
477  niter = 1
478  ip1 = i + 1
479 *
480 * Calculate initial guess
481 *
482  del = d( ip1 ) - d( i )
483  midpt = del / two
484  DO 100 j = 1, n
485  delta( j ) = ( d( j )-d( i ) ) - midpt
486  100 CONTINUE
487 *
488  psi = zero
489  DO 110 j = 1, i - 1
490  psi = psi + z( j )*z( j ) / delta( j )
491  110 CONTINUE
492 *
493  phi = zero
494  DO 120 j = n, i + 2, -1
495  phi = phi + z( j )*z( j ) / delta( j )
496  120 CONTINUE
497  c = rhoinv + psi + phi
498  w = c + z( i )*z( i ) / delta( i ) +
499  $ z( ip1 )*z( ip1 ) / delta( ip1 )
500 *
501  IF( w.GT.zero ) THEN
502 *
503 * d(i)< the ith eigenvalue < (d(i)+d(i+1))/2
504 *
505 * We choose d(i) as origin.
506 *
507  orgati = .true.
508  a = c*del + z( i )*z( i ) + z( ip1 )*z( ip1 )
509  b = z( i )*z( i )*del
510  IF( a.GT.zero ) THEN
511  tau = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
512  ELSE
513  tau = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
514  END IF
515  dltlb = zero
516  dltub = midpt
517  ELSE
518 *
519 * (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1)
520 *
521 * We choose d(i+1) as origin.
522 *
523  orgati = .false.
524  a = c*del - z( i )*z( i ) - z( ip1 )*z( ip1 )
525  b = z( ip1 )*z( ip1 )*del
526  IF( a.LT.zero ) THEN
527  tau = two*b / ( a-sqrt( abs( a*a+four*b*c ) ) )
528  ELSE
529  tau = -( a+sqrt( abs( a*a+four*b*c ) ) ) / ( two*c )
530  END IF
531  dltlb = -midpt
532  dltub = zero
533  END IF
534 *
535  IF( orgati ) THEN
536  DO 130 j = 1, n
537  delta( j ) = ( d( j )-d( i ) ) - tau
538  130 CONTINUE
539  ELSE
540  DO 140 j = 1, n
541  delta( j ) = ( d( j )-d( ip1 ) ) - tau
542  140 CONTINUE
543  END IF
544  IF( orgati ) THEN
545  ii = i
546  ELSE
547  ii = i + 1
548  END IF
549  iim1 = ii - 1
550  iip1 = ii + 1
551 *
552 * Evaluate PSI and the derivative DPSI
553 *
554  dpsi = zero
555  psi = zero
556  erretm = zero
557  DO 150 j = 1, iim1
558  temp = z( j ) / delta( j )
559  psi = psi + z( j )*temp
560  dpsi = dpsi + temp*temp
561  erretm = erretm + psi
562  150 CONTINUE
563  erretm = abs( erretm )
564 *
565 * Evaluate PHI and the derivative DPHI
566 *
567  dphi = zero
568  phi = zero
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
574  160 CONTINUE
575 *
576  w = rhoinv + phi + psi
577 *
578 * W is the value of the secular function with
579 * its ii-th element removed.
580 *
581  swtch3 = .false.
582  IF( orgati ) THEN
583  IF( w.LT.zero )
584  $ swtch3 = .true.
585  ELSE
586  IF( w.GT.zero )
587  $ swtch3 = .true.
588  END IF
589  IF( ii.EQ.1 .OR. ii.EQ.n )
590  $ swtch3 = .false.
591 *
592  temp = z( ii ) / delta( ii )
593  dw = dpsi + dphi + temp*temp
594  temp = z( ii )*temp
595  w = w + temp
596  erretm = eight*( phi-psi ) + erretm + two*rhoinv +
597  $ three*abs( temp ) + abs( tau )*dw
598 *
599 * Test for convergence
600 *
601  IF( abs( w ).LE.eps*erretm ) THEN
602  IF( orgati ) THEN
603  dlam = d( i ) + tau
604  ELSE
605  dlam = d( ip1 ) + tau
606  END IF
607  GO TO 250
608  END IF
609 *
610  IF( w.LE.zero ) THEN
611  dltlb = max( dltlb, tau )
612  ELSE
613  dltub = min( dltub, tau )
614  END IF
615 *
616 * Calculate the new step
617 *
618  niter = niter + 1
619  IF( .NOT.swtch3 ) THEN
620  IF( orgati ) THEN
621  c = w - delta( ip1 )*dw - ( d( i )-d( ip1 ) )*
622  $ ( z( i ) / delta( i ) )**2
623  ELSE
624  c = w - delta( i )*dw - ( d( ip1 )-d( i ) )*
625  $ ( z( ip1 ) / delta( ip1 ) )**2
626  END IF
627  a = ( delta( i )+delta( ip1 ) )*w -
628  $ delta( i )*delta( ip1 )*dw
629  b = delta( i )*delta( ip1 )*w
630  IF( c.EQ.zero ) THEN
631  IF( a.EQ.zero ) THEN
632  IF( orgati ) THEN
633  a = z( i )*z( i ) + delta( ip1 )*delta( ip1 )*
634  $ ( dpsi+dphi )
635  ELSE
636  a = z( ip1 )*z( ip1 ) + delta( i )*delta( i )*
637  $ ( dpsi+dphi )
638  END IF
639  END IF
640  eta = b / a
641  ELSE IF( a.LE.zero ) THEN
642  eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
643  ELSE
644  eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
645  END IF
646  ELSE
647 *
648 * Interpolation using THREE most relevant poles
649 *
650  temp = rhoinv + psi + phi
651  IF( orgati ) THEN
652  temp1 = z( iim1 ) / delta( iim1 )
653  temp1 = temp1*temp1
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 )
659  ELSE
660  temp1 = z( iip1 ) / delta( iip1 )
661  temp1 = temp1*temp1
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 )
667  END IF
668  zz( 2 ) = z( ii )*z( ii )
669  CALL slaed6( niter, orgati, c, delta( iim1 ), zz, w, eta,
670  $ info )
671  IF( info.NE.0 )
672  $ GO TO 250
673  END IF
674 *
675 * Note, eta should be positive if w is negative, and
676 * eta should be negative otherwise. However,
677 * if for some reason caused by roundoff, eta*w > 0,
678 * we simply use one Newton step instead. This way
679 * will guarantee eta*w < 0.
680 *
681  IF( w*eta.GE.zero )
682  $ eta = -w / dw
683  temp = tau + eta
684  IF( temp.GT.dltub .OR. temp.LT.dltlb ) THEN
685  IF( w.LT.zero ) THEN
686  eta = ( dltub-tau ) / two
687  ELSE
688  eta = ( dltlb-tau ) / two
689  END IF
690  END IF
691 *
692  prew = w
693 *
694  DO 180 j = 1, n
695  delta( j ) = delta( j ) - eta
696  180 CONTINUE
697 *
698 * Evaluate PSI and the derivative DPSI
699 *
700  dpsi = zero
701  psi = zero
702  erretm = zero
703  DO 190 j = 1, iim1
704  temp = z( j ) / delta( j )
705  psi = psi + z( j )*temp
706  dpsi = dpsi + temp*temp
707  erretm = erretm + psi
708  190 CONTINUE
709  erretm = abs( erretm )
710 *
711 * Evaluate PHI and the derivative DPHI
712 *
713  dphi = zero
714  phi = zero
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
720  200 CONTINUE
721 *
722  temp = z( ii ) / delta( ii )
723  dw = dpsi + dphi + temp*temp
724  temp = z( ii )*temp
725  w = rhoinv + phi + psi + temp
726  erretm = eight*( phi-psi ) + erretm + two*rhoinv +
727  $ three*abs( temp ) + abs( tau+eta )*dw
728 *
729  swtch = .false.
730  IF( orgati ) THEN
731  IF( -w.GT.abs( prew ) / ten )
732  $ swtch = .true.
733  ELSE
734  IF( w.GT.abs( prew ) / ten )
735  $ swtch = .true.
736  END IF
737 *
738  tau = tau + eta
739 *
740 * Main loop to update the values of the array DELTA
741 *
742  iter = niter + 1
743 *
744  DO 240 niter = iter, maxit
745 *
746 * Test for convergence
747 *
748  IF( abs( w ).LE.eps*erretm ) THEN
749  IF( orgati ) THEN
750  dlam = d( i ) + tau
751  ELSE
752  dlam = d( ip1 ) + tau
753  END IF
754  GO TO 250
755  END IF
756 *
757  IF( w.LE.zero ) THEN
758  dltlb = max( dltlb, tau )
759  ELSE
760  dltub = min( dltub, tau )
761  END IF
762 *
763 * Calculate the new step
764 *
765  IF( .NOT.swtch3 ) THEN
766  IF( .NOT.swtch ) THEN
767  IF( orgati ) THEN
768  c = w - delta( ip1 )*dw -
769  $ ( d( i )-d( ip1 ) )*( z( i ) / delta( i ) )**2
770  ELSE
771  c = w - delta( i )*dw - ( d( ip1 )-d( i ) )*
772  $ ( z( ip1 ) / delta( ip1 ) )**2
773  END IF
774  ELSE
775  temp = z( ii ) / delta( ii )
776  IF( orgati ) THEN
777  dpsi = dpsi + temp*temp
778  ELSE
779  dphi = dphi + temp*temp
780  END IF
781  c = w - delta( i )*dpsi - delta( ip1 )*dphi
782  END IF
783  a = ( delta( i )+delta( ip1 ) )*w -
784  $ delta( i )*delta( ip1 )*dw
785  b = delta( i )*delta( ip1 )*w
786  IF( c.EQ.zero ) THEN
787  IF( a.EQ.zero ) THEN
788  IF( .NOT.swtch ) THEN
789  IF( orgati ) THEN
790  a = z( i )*z( i ) + delta( ip1 )*
791  $ delta( ip1 )*( dpsi+dphi )
792  ELSE
793  a = z( ip1 )*z( ip1 ) +
794  $ delta( i )*delta( i )*( dpsi+dphi )
795  END IF
796  ELSE
797  a = delta( i )*delta( i )*dpsi +
798  $ delta( ip1 )*delta( ip1 )*dphi
799  END IF
800  END IF
801  eta = b / a
802  ELSE IF( a.LE.zero ) THEN
803  eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
804  ELSE
805  eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
806  END IF
807  ELSE
808 *
809 * Interpolation using THREE most relevant poles
810 *
811  temp = rhoinv + psi + phi
812  IF( swtch ) THEN
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
816  ELSE
817  IF( orgati ) THEN
818  temp1 = z( iim1 ) / delta( iim1 )
819  temp1 = temp1*temp1
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 )
825  ELSE
826  temp1 = z( iip1 ) / delta( iip1 )
827  temp1 = temp1*temp1
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 )
833  END IF
834  END IF
835  CALL slaed6( niter, orgati, c, delta( iim1 ), zz, w, eta,
836  $ info )
837  IF( info.NE.0 )
838  $ GO TO 250
839  END IF
840 *
841 * Note, eta should be positive if w is negative, and
842 * eta should be negative otherwise. However,
843 * if for some reason caused by roundoff, eta*w > 0,
844 * we simply use one Newton step instead. This way
845 * will guarantee eta*w < 0.
846 *
847  IF( w*eta.GE.zero )
848  $ eta = -w / dw
849  temp = tau + eta
850  IF( temp.GT.dltub .OR. temp.LT.dltlb ) THEN
851  IF( w.LT.zero ) THEN
852  eta = ( dltub-tau ) / two
853  ELSE
854  eta = ( dltlb-tau ) / two
855  END IF
856  END IF
857 *
858  DO 210 j = 1, n
859  delta( j ) = delta( j ) - eta
860  210 CONTINUE
861 *
862  tau = tau + eta
863  prew = w
864 *
865 * Evaluate PSI and the derivative DPSI
866 *
867  dpsi = zero
868  psi = zero
869  erretm = zero
870  DO 220 j = 1, iim1
871  temp = z( j ) / delta( j )
872  psi = psi + z( j )*temp
873  dpsi = dpsi + temp*temp
874  erretm = erretm + psi
875  220 CONTINUE
876  erretm = abs( erretm )
877 *
878 * Evaluate PHI and the derivative DPHI
879 *
880  dphi = zero
881  phi = zero
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
887  230 CONTINUE
888 *
889  temp = z( ii ) / delta( ii )
890  dw = dpsi + dphi + temp*temp
891  temp = z( ii )*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 )
896  $ swtch = .NOT.swtch
897 *
898  240 CONTINUE
899 *
900 * Return with INFO = 1, NITER = MAXIT and not converged
901 *
902  info = 1
903  IF( orgati ) THEN
904  dlam = d( i ) + tau
905  ELSE
906  dlam = d( ip1 ) + tau
907  END IF
908 *
909  END IF
910 *
911  250 CONTINUE
912 *
913  RETURN
914 *
915 * End of SLAED4
916 *
slaed5
subroutine slaed5(I, D, Z, DELTA, RHO, DLAM)
SLAED5 used by sstedc. Solves the 2-by-2 secular equation.
Definition: slaed5.f:110
slaed6
subroutine slaed6(KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)
SLAED6 used by sstedc. Computes one Newton step in solution of the secular equation.
Definition: slaed6.f:142
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69

Here is the call graph for this function:

Here is the caller graph for this function: