LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlaed4()

subroutine dlaed4 ( integer  n,
integer  i,
double precision, dimension( * )  d,
double precision, dimension( * )  z,
double precision, dimension( * )  delta,
double precision  rho,
double precision  dlam,
integer  info 
)

DLAED4 used by DSTEDC. Finds a single root of the secular equation.

Download DLAED4 + 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
         The components of the updating vector.
[out]DELTA
          DELTA is DOUBLE PRECISION 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 DLAED5
         for detail. The vector DELTA contains the information necessary
         to construct the eigenvectors by DLAED3 and DLAED9.
[in]RHO
          RHO is DOUBLE PRECISION
         The scalar in the symmetric updating formula.
[out]DLAM
          DLAM is DOUBLE PRECISION
         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 dlaed4.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 DOUBLE PRECISION DLAM, RHO
153* ..
154* .. Array Arguments ..
155 DOUBLE PRECISION D( * ), DELTA( * ), Z( * )
156* ..
157*
158* =====================================================================
159*
160* .. Parameters ..
161 INTEGER MAXIT
162 parameter( maxit = 30 )
163 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
164 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
165 $ three = 3.0d0, four = 4.0d0, eight = 8.0d0,
166 $ ten = 10.0d0 )
167* ..
168* .. Local Scalars ..
169 LOGICAL ORGATI, SWTCH, SWTCH3
170 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
171 DOUBLE PRECISION 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 DOUBLE PRECISION ZZ( 3 )
177* ..
178* .. External Functions ..
179 DOUBLE PRECISION DLAMCH
180 EXTERNAL dlamch
181* ..
182* .. External Subroutines ..
183 EXTERNAL dlaed5, dlaed6
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 dlaed5( i, d, z, delta, rho, dlam )
206 RETURN
207 END IF
208*
209* Compute machine epsilon
210*
211 eps = dlamch( '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*
335* Update proposed by Li, Ren-Cang:
336 eta = -w / ( dpsi+dphi )
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 dlaed6( 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 dlaed6( 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 DLAED4
916*
dlaed5
subroutine dlaed5(i, d, z, delta, rho, dlam)
DLAED5 used by DSTEDC. Solves the 2-by-2 secular equation.
Definition dlaed5.f:108
dlaed6
subroutine dlaed6(kniter, orgati, rho, d, z, finit, tau, info)
DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation.
Definition dlaed6.f:140
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
Here is the call graph for this function:
Here is the caller graph for this function: