LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slaed4()

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 > 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.
Contributors:
Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

Definition at line 144 of file slaed4.f.

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