LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaed4.f
Go to the documentation of this file.
1*> \brief \b DLAED4 used by DSTEDC. Finds a single root of the secular equation.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DLAED4 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed4.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed4.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed4.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
20*
21* .. Scalar Arguments ..
22* INTEGER I, INFO, N
23* DOUBLE PRECISION DLAM, RHO
24* ..
25* .. Array Arguments ..
26* DOUBLE PRECISION D( * ), DELTA( * ), Z( * )
27* ..
28*
29*
30*> \par Purpose:
31* =============
32*>
33*> \verbatim
34*>
35*> This subroutine computes the I-th updated eigenvalue of a symmetric
36*> rank-one modification to a diagonal matrix whose elements are
37*> given in the array d, and that
38*>
39*> D(i) < D(j) for i < j
40*>
41*> and that RHO > 0. This is arranged by the calling routine, and is
42*> no loss in generality. The rank-one modified system is thus
43*>
44*> diag( D ) + RHO * Z * Z_transpose.
45*>
46*> where we assume the Euclidean norm of Z is 1.
47*>
48*> The method consists of approximating the rational functions in the
49*> secular equation by simpler interpolating rational functions.
50*> \endverbatim
51*
52* Arguments:
53* ==========
54*
55*> \param[in] N
56*> \verbatim
57*> N is INTEGER
58*> The length of all arrays.
59*> \endverbatim
60*>
61*> \param[in] I
62*> \verbatim
63*> I is INTEGER
64*> The index of the eigenvalue to be computed. 1 <= I <= N.
65*> \endverbatim
66*>
67*> \param[in] D
68*> \verbatim
69*> D is DOUBLE PRECISION array, dimension (N)
70*> The original eigenvalues. It is assumed that they are in
71*> order, D(I) < D(J) for I < J.
72*> \endverbatim
73*>
74*> \param[in] Z
75*> \verbatim
76*> Z is DOUBLE PRECISION array, dimension (N)
77*> The components of the updating vector.
78*> \endverbatim
79*>
80*> \param[out] DELTA
81*> \verbatim
82*> DELTA is DOUBLE PRECISION array, dimension (N)
83*> If N > 2, DELTA contains (D(j) - lambda_I) in its j-th
84*> component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5
85*> for detail. The vector DELTA contains the information necessary
86*> to construct the eigenvectors by DLAED3 and DLAED9.
87*> \endverbatim
88*>
89*> \param[in] RHO
90*> \verbatim
91*> RHO is DOUBLE PRECISION
92*> The scalar in the symmetric updating formula.
93*> \endverbatim
94*>
95*> \param[out] DLAM
96*> \verbatim
97*> DLAM is DOUBLE PRECISION
98*> The computed lambda_I, the I-th updated eigenvalue.
99*> \endverbatim
100*>
101*> \param[out] INFO
102*> \verbatim
103*> INFO is INTEGER
104*> = 0: successful exit
105*> > 0: if INFO = 1, the updating process failed.
106*> \endverbatim
107*
108*> \par Internal Parameters:
109* =========================
110*>
111*> \verbatim
112*> Logical variable ORGATI (origin-at-i?) is used for distinguishing
113*> whether D(i) or D(i+1) is treated as the origin.
114*>
115*> ORGATI = .true. origin at i
116*> ORGATI = .false. origin at i+1
117*>
118*> Logical variable SWTCH3 (switch-for-3-poles?) is for noting
119*> if we are working with THREE poles!
120*>
121*> MAXIT is the maximum number of iterations allowed for each
122*> eigenvalue.
123*> \endverbatim
124*
125* Authors:
126* ========
127*
128*> \author Univ. of Tennessee
129*> \author Univ. of California Berkeley
130*> \author Univ. of Colorado Denver
131*> \author NAG Ltd.
132*
133*> \ingroup laed4
134*
135*> \par Contributors:
136* ==================
137*>
138*> Ren-Cang Li, Computer Science Division, University of California
139*> at Berkeley, USA
140*>
141* =====================================================================
142 SUBROUTINE dlaed4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
143*
144* -- LAPACK computational routine --
145* -- LAPACK is a software package provided by Univ. of Tennessee, --
146* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
147*
148* .. Scalar Arguments ..
149 INTEGER I, INFO, N
150 DOUBLE PRECISION DLAM, RHO
151* ..
152* .. Array Arguments ..
153 DOUBLE PRECISION D( * ), DELTA( * ), Z( * )
154* ..
155*
156* =====================================================================
157*
158* .. Parameters ..
159 INTEGER MAXIT
160 parameter( maxit = 30 )
161 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
162 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
163 $ three = 3.0d0, four = 4.0d0, eight = 8.0d0,
164 $ ten = 10.0d0 )
165* ..
166* .. Local Scalars ..
167 LOGICAL ORGATI, SWTCH, SWTCH3
168 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
169 DOUBLE PRECISION A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
170 $ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
171 $ RHOINV, TAU, TEMP, TEMP1, W
172* ..
173* .. Local Arrays ..
174 DOUBLE PRECISION ZZ( 3 )
175* ..
176* .. External Functions ..
177 DOUBLE PRECISION DLAMCH
178 EXTERNAL dlamch
179* ..
180* .. External Subroutines ..
181 EXTERNAL dlaed5, dlaed6
182* ..
183* .. Intrinsic Functions ..
184 INTRINSIC abs, max, min, sqrt
185* ..
186* .. Executable Statements ..
187*
188* Since this routine is called in an inner loop, we do no argument
189* checking.
190*
191* Quick return for N=1 and 2.
192*
193 info = 0
194 IF( n.EQ.1 ) THEN
195*
196* Presumably, I=1 upon entry
197*
198 dlam = d( 1 ) + rho*z( 1 )*z( 1 )
199 delta( 1 ) = one
200 RETURN
201 END IF
202 IF( n.EQ.2 ) THEN
203 CALL dlaed5( i, d, z, delta, rho, dlam )
204 RETURN
205 END IF
206*
207* Compute machine epsilon
208*
209 eps = dlamch( 'Epsilon' )
210 rhoinv = one / rho
211*
212* The case I = N
213*
214 IF( i.EQ.n ) THEN
215*
216* Initialize some basic variables
217*
218 ii = n - 1
219 niter = 1
220*
221* Calculate initial guess
222*
223 midpt = rho / two
224*
225* If ||Z||_2 is not one, then TEMP should be set to
226* RHO * ||Z||_2^2 / TWO
227*
228 DO 10 j = 1, n
229 delta( j ) = ( d( j )-d( i ) ) - midpt
230 10 CONTINUE
231*
232 psi = zero
233 DO 20 j = 1, n - 2
234 psi = psi + z( j )*z( j ) / delta( j )
235 20 CONTINUE
236*
237 c = rhoinv + psi
238 w = c + z( ii )*z( ii ) / delta( ii ) +
239 $ z( n )*z( n ) / delta( n )
240*
241 IF( w.LE.zero ) THEN
242 temp = z( n-1 )*z( n-1 ) / ( d( n )-d( n-1 )+rho ) +
243 $ z( n )*z( n ) / rho
244 IF( c.LE.temp ) THEN
245 tau = rho
246 ELSE
247 del = d( n ) - d( n-1 )
248 a = -c*del + z( n-1 )*z( n-1 ) + z( n )*z( n )
249 b = z( n )*z( n )*del
250 IF( a.LT.zero ) THEN
251 tau = two*b / ( sqrt( a*a+four*b*c )-a )
252 ELSE
253 tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
254 END IF
255 END IF
256*
257* It can be proved that
258* D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO
259*
260 dltlb = midpt
261 dltub = rho
262 ELSE
263 del = d( n ) - d( n-1 )
264 a = -c*del + z( n-1 )*z( n-1 ) + z( n )*z( n )
265 b = z( n )*z( n )*del
266 IF( a.LT.zero ) THEN
267 tau = two*b / ( sqrt( a*a+four*b*c )-a )
268 ELSE
269 tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
270 END IF
271*
272* It can be proved that
273* D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2
274*
275 dltlb = zero
276 dltub = midpt
277 END IF
278*
279 DO 30 j = 1, n
280 delta( j ) = ( d( j )-d( i ) ) - tau
281 30 CONTINUE
282*
283* Evaluate PSI and the derivative DPSI
284*
285 dpsi = zero
286 psi = zero
287 erretm = zero
288 DO 40 j = 1, ii
289 temp = z( j ) / delta( j )
290 psi = psi + z( j )*temp
291 dpsi = dpsi + temp*temp
292 erretm = erretm + psi
293 40 CONTINUE
294 erretm = abs( erretm )
295*
296* Evaluate PHI and the derivative DPHI
297*
298 temp = z( n ) / delta( n )
299 phi = z( n )*temp
300 dphi = temp*temp
301 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
302 $ abs( tau )*( dpsi+dphi )
303*
304 w = rhoinv + phi + psi
305*
306* Test for convergence
307*
308 IF( abs( w ).LE.eps*erretm ) THEN
309 dlam = d( i ) + tau
310 GO TO 250
311 END IF
312*
313 IF( w.LE.zero ) THEN
314 dltlb = max( dltlb, tau )
315 ELSE
316 dltub = min( dltub, tau )
317 END IF
318*
319* Calculate the new step
320*
321 niter = niter + 1
322 c = w - delta( n-1 )*dpsi - delta( n )*dphi
323 a = ( delta( n-1 )+delta( n ) )*w -
324 $ delta( n-1 )*delta( n )*( dpsi+dphi )
325 b = delta( n-1 )*delta( n )*w
326 IF( c.LT.zero )
327 $ c = abs( c )
328 IF( c.EQ.zero ) THEN
329* ETA = B/A
330* ETA = RHO - TAU
331* ETA = DLTUB - TAU
332*
333* Update proposed by Li, Ren-Cang:
334 eta = -w / ( dpsi+dphi )
335 ELSE IF( a.GE.zero ) THEN
336 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
337 ELSE
338 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
339 END IF
340*
341* Note, eta should be positive if w is negative, and
342* eta should be negative otherwise. However,
343* if for some reason caused by roundoff, eta*w > 0,
344* we simply use one Newton step instead. This way
345* will guarantee eta*w < 0.
346*
347 IF( w*eta.GT.zero )
348 $ eta = -w / ( dpsi+dphi )
349 temp = tau + eta
350 IF( temp.GT.dltub .OR. temp.LT.dltlb ) THEN
351 IF( w.LT.zero ) THEN
352 eta = ( dltub-tau ) / two
353 ELSE
354 eta = ( dltlb-tau ) / two
355 END IF
356 END IF
357 DO 50 j = 1, n
358 delta( j ) = delta( j ) - eta
359 50 CONTINUE
360*
361 tau = tau + eta
362*
363* Evaluate PSI and the derivative DPSI
364*
365 dpsi = zero
366 psi = zero
367 erretm = zero
368 DO 60 j = 1, ii
369 temp = z( j ) / delta( j )
370 psi = psi + z( j )*temp
371 dpsi = dpsi + temp*temp
372 erretm = erretm + psi
373 60 CONTINUE
374 erretm = abs( erretm )
375*
376* Evaluate PHI and the derivative DPHI
377*
378 temp = z( n ) / delta( n )
379 phi = z( n )*temp
380 dphi = temp*temp
381 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
382 $ abs( tau )*( dpsi+dphi )
383*
384 w = rhoinv + phi + psi
385*
386* Main loop to update the values of the array DELTA
387*
388 iter = niter + 1
389*
390 DO 90 niter = iter, maxit
391*
392* Test for convergence
393*
394 IF( abs( w ).LE.eps*erretm ) THEN
395 dlam = d( i ) + tau
396 GO TO 250
397 END IF
398*
399 IF( w.LE.zero ) THEN
400 dltlb = max( dltlb, tau )
401 ELSE
402 dltub = min( dltub, tau )
403 END IF
404*
405* Calculate the new step
406*
407 c = w - delta( n-1 )*dpsi - delta( n )*dphi
408 a = ( delta( n-1 )+delta( n ) )*w -
409 $ delta( n-1 )*delta( n )*( dpsi+dphi )
410 b = delta( n-1 )*delta( n )*w
411 IF( a.GE.zero ) THEN
412 eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
413 ELSE
414 eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
415 END IF
416*
417* Note, eta should be positive if w is negative, and
418* eta should be negative otherwise. However,
419* if for some reason caused by roundoff, eta*w > 0,
420* we simply use one Newton step instead. This way
421* will guarantee eta*w < 0.
422*
423 IF( w*eta.GT.zero )
424 $ eta = -w / ( dpsi+dphi )
425 temp = tau + eta
426 IF( temp.GT.dltub .OR. temp.LT.dltlb ) THEN
427 IF( w.LT.zero ) THEN
428 eta = ( dltub-tau ) / two
429 ELSE
430 eta = ( dltlb-tau ) / two
431 END IF
432 END IF
433 DO 70 j = 1, n
434 delta( j ) = delta( j ) - eta
435 70 CONTINUE
436*
437 tau = tau + eta
438*
439* Evaluate PSI and the derivative DPSI
440*
441 dpsi = zero
442 psi = zero
443 erretm = zero
444 DO 80 j = 1, ii
445 temp = z( j ) / delta( j )
446 psi = psi + z( j )*temp
447 dpsi = dpsi + temp*temp
448 erretm = erretm + psi
449 80 CONTINUE
450 erretm = abs( erretm )
451*
452* Evaluate PHI and the derivative DPHI
453*
454 temp = z( n ) / delta( n )
455 phi = z( n )*temp
456 dphi = temp*temp
457 erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
458 $ abs( tau )*( dpsi+dphi )
459*
460 w = rhoinv + phi + psi
461 90 CONTINUE
462*
463* Return with INFO = 1, NITER = MAXIT and not converged
464*
465 info = 1
466 dlam = d( i ) + tau
467 GO TO 250
468*
469* End for the case I = N
470*
471 ELSE
472*
473* The case for I < N
474*
475 niter = 1
476 ip1 = i + 1
477*
478* Calculate initial guess
479*
480 del = d( ip1 ) - d( i )
481 midpt = del / two
482 DO 100 j = 1, n
483 delta( j ) = ( d( j )-d( i ) ) - midpt
484 100 CONTINUE
485*
486 psi = zero
487 DO 110 j = 1, i - 1
488 psi = psi + z( j )*z( j ) / delta( j )
489 110 CONTINUE
490*
491 phi = zero
492 DO 120 j = n, i + 2, -1
493 phi = phi + z( j )*z( j ) / delta( j )
494 120 CONTINUE
495 c = rhoinv + psi + phi
496 w = c + z( i )*z( i ) / delta( i ) +
497 $ z( ip1 )*z( ip1 ) / delta( ip1 )
498*
499 IF( w.GT.zero ) THEN
500*
501* d(i)< the ith eigenvalue < (d(i)+d(i+1))/2
502*
503* We choose d(i) as origin.
504*
505 orgati = .true.
506 a = c*del + z( i )*z( i ) + z( ip1 )*z( ip1 )
507 b = z( i )*z( i )*del
508 IF( a.GT.zero ) THEN
509 tau = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
510 ELSE
511 tau = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
512 END IF
513 dltlb = zero
514 dltub = midpt
515 ELSE
516*
517* (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1)
518*
519* We choose d(i+1) as origin.
520*
521 orgati = .false.
522 a = c*del - z( i )*z( i ) - z( ip1 )*z( ip1 )
523 b = z( ip1 )*z( ip1 )*del
524 IF( a.LT.zero ) THEN
525 tau = two*b / ( a-sqrt( abs( a*a+four*b*c ) ) )
526 ELSE
527 tau = -( a+sqrt( abs( a*a+four*b*c ) ) ) / ( two*c )
528 END IF
529 dltlb = -midpt
530 dltub = zero
531 END IF
532*
533 IF( orgati ) THEN
534 DO 130 j = 1, n
535 delta( j ) = ( d( j )-d( i ) ) - tau
536 130 CONTINUE
537 ELSE
538 DO 140 j = 1, n
539 delta( j ) = ( d( j )-d( ip1 ) ) - tau
540 140 CONTINUE
541 END IF
542 IF( orgati ) THEN
543 ii = i
544 ELSE
545 ii = i + 1
546 END IF
547 iim1 = ii - 1
548 iip1 = ii + 1
549*
550* Evaluate PSI and the derivative DPSI
551*
552 dpsi = zero
553 psi = zero
554 erretm = zero
555 DO 150 j = 1, iim1
556 temp = z( j ) / delta( j )
557 psi = psi + z( j )*temp
558 dpsi = dpsi + temp*temp
559 erretm = erretm + psi
560 150 CONTINUE
561 erretm = abs( erretm )
562*
563* Evaluate PHI and the derivative DPHI
564*
565 dphi = zero
566 phi = zero
567 DO 160 j = n, iip1, -1
568 temp = z( j ) / delta( j )
569 phi = phi + z( j )*temp
570 dphi = dphi + temp*temp
571 erretm = erretm + phi
572 160 CONTINUE
573*
574 w = rhoinv + phi + psi
575*
576* W is the value of the secular function with
577* its ii-th element removed.
578*
579 swtch3 = .false.
580 IF( orgati ) THEN
581 IF( w.LT.zero )
582 $ swtch3 = .true.
583 ELSE
584 IF( w.GT.zero )
585 $ swtch3 = .true.
586 END IF
587 IF( ii.EQ.1 .OR. ii.EQ.n )
588 $ swtch3 = .false.
589*
590 temp = z( ii ) / delta( ii )
591 dw = dpsi + dphi + temp*temp
592 temp = z( ii )*temp
593 w = w + temp
594 erretm = eight*( phi-psi ) + erretm + two*rhoinv +
595 $ three*abs( temp ) + abs( tau )*dw
596*
597* Test for convergence
598*
599 IF( abs( w ).LE.eps*erretm ) THEN
600 IF( orgati ) THEN
601 dlam = d( i ) + tau
602 ELSE
603 dlam = d( ip1 ) + tau
604 END IF
605 GO TO 250
606 END IF
607*
608 IF( w.LE.zero ) THEN
609 dltlb = max( dltlb, tau )
610 ELSE
611 dltub = min( dltub, tau )
612 END IF
613*
614* Calculate the new step
615*
616 niter = niter + 1
617 IF( .NOT.swtch3 ) THEN
618 IF( orgati ) THEN
619 c = w - delta( ip1 )*dw - ( d( i )-d( ip1 ) )*
620 $ ( z( i ) / delta( i ) )**2
621 ELSE
622 c = w - delta( i )*dw - ( d( ip1 )-d( i ) )*
623 $ ( z( ip1 ) / delta( ip1 ) )**2
624 END IF
625 a = ( delta( i )+delta( ip1 ) )*w -
626 $ delta( i )*delta( ip1 )*dw
627 b = delta( i )*delta( ip1 )*w
628 IF( c.EQ.zero ) THEN
629 IF( a.EQ.zero ) THEN
630 IF( orgati ) THEN
631 a = z( i )*z( i ) + delta( ip1 )*delta( ip1 )*
632 $ ( dpsi+dphi )
633 ELSE
634 a = z( ip1 )*z( ip1 ) + delta( i )*delta( i )*
635 $ ( dpsi+dphi )
636 END IF
637 END IF
638 eta = b / a
639 ELSE IF( a.LE.zero ) THEN
640 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
641 ELSE
642 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
643 END IF
644 ELSE
645*
646* Interpolation using THREE most relevant poles
647*
648 temp = rhoinv + psi + phi
649 IF( orgati ) THEN
650 temp1 = z( iim1 ) / delta( iim1 )
651 temp1 = temp1*temp1
652 c = temp - delta( iip1 )*( dpsi+dphi ) -
653 $ ( d( iim1 )-d( iip1 ) )*temp1
654 zz( 1 ) = z( iim1 )*z( iim1 )
655 zz( 3 ) = delta( iip1 )*delta( iip1 )*
656 $ ( ( dpsi-temp1 )+dphi )
657 ELSE
658 temp1 = z( iip1 ) / delta( iip1 )
659 temp1 = temp1*temp1
660 c = temp - delta( iim1 )*( dpsi+dphi ) -
661 $ ( d( iip1 )-d( iim1 ) )*temp1
662 zz( 1 ) = delta( iim1 )*delta( iim1 )*
663 $ ( dpsi+( dphi-temp1 ) )
664 zz( 3 ) = z( iip1 )*z( iip1 )
665 END IF
666 zz( 2 ) = z( ii )*z( ii )
667 CALL dlaed6( niter, orgati, c, delta( iim1 ), zz, w, eta,
668 $ info )
669 IF( info.NE.0 )
670 $ GO TO 250
671 END IF
672*
673* Note, eta should be positive if w is negative, and
674* eta should be negative otherwise. However,
675* if for some reason caused by roundoff, eta*w > 0,
676* we simply use one Newton step instead. This way
677* will guarantee eta*w < 0.
678*
679 IF( w*eta.GE.zero )
680 $ eta = -w / dw
681 temp = tau + eta
682 IF( temp.GT.dltub .OR. temp.LT.dltlb ) THEN
683 IF( w.LT.zero ) THEN
684 eta = ( dltub-tau ) / two
685 ELSE
686 eta = ( dltlb-tau ) / two
687 END IF
688 END IF
689*
690 prew = w
691*
692 DO 180 j = 1, n
693 delta( j ) = delta( j ) - eta
694 180 CONTINUE
695*
696* Evaluate PSI and the derivative DPSI
697*
698 dpsi = zero
699 psi = zero
700 erretm = zero
701 DO 190 j = 1, iim1
702 temp = z( j ) / delta( j )
703 psi = psi + z( j )*temp
704 dpsi = dpsi + temp*temp
705 erretm = erretm + psi
706 190 CONTINUE
707 erretm = abs( erretm )
708*
709* Evaluate PHI and the derivative DPHI
710*
711 dphi = zero
712 phi = zero
713 DO 200 j = n, iip1, -1
714 temp = z( j ) / delta( j )
715 phi = phi + z( j )*temp
716 dphi = dphi + temp*temp
717 erretm = erretm + phi
718 200 CONTINUE
719*
720 temp = z( ii ) / delta( ii )
721 dw = dpsi + dphi + temp*temp
722 temp = z( ii )*temp
723 w = rhoinv + phi + psi + temp
724 erretm = eight*( phi-psi ) + erretm + two*rhoinv +
725 $ three*abs( temp ) + abs( tau+eta )*dw
726*
727 swtch = .false.
728 IF( orgati ) THEN
729 IF( -w.GT.abs( prew ) / ten )
730 $ swtch = .true.
731 ELSE
732 IF( w.GT.abs( prew ) / ten )
733 $ swtch = .true.
734 END IF
735*
736 tau = tau + eta
737*
738* Main loop to update the values of the array DELTA
739*
740 iter = niter + 1
741*
742 DO 240 niter = iter, maxit
743*
744* Test for convergence
745*
746 IF( abs( w ).LE.eps*erretm ) THEN
747 IF( orgati ) THEN
748 dlam = d( i ) + tau
749 ELSE
750 dlam = d( ip1 ) + tau
751 END IF
752 GO TO 250
753 END IF
754*
755 IF( w.LE.zero ) THEN
756 dltlb = max( dltlb, tau )
757 ELSE
758 dltub = min( dltub, tau )
759 END IF
760*
761* Calculate the new step
762*
763 IF( .NOT.swtch3 ) THEN
764 IF( .NOT.swtch ) THEN
765 IF( orgati ) THEN
766 c = w - delta( ip1 )*dw -
767 $ ( d( i )-d( ip1 ) )*( z( i ) / delta( i ) )**2
768 ELSE
769 c = w - delta( i )*dw - ( d( ip1 )-d( i ) )*
770 $ ( z( ip1 ) / delta( ip1 ) )**2
771 END IF
772 ELSE
773 temp = z( ii ) / delta( ii )
774 IF( orgati ) THEN
775 dpsi = dpsi + temp*temp
776 ELSE
777 dphi = dphi + temp*temp
778 END IF
779 c = w - delta( i )*dpsi - delta( ip1 )*dphi
780 END IF
781 a = ( delta( i )+delta( ip1 ) )*w -
782 $ delta( i )*delta( ip1 )*dw
783 b = delta( i )*delta( ip1 )*w
784 IF( c.EQ.zero ) THEN
785 IF( a.EQ.zero ) THEN
786 IF( .NOT.swtch ) THEN
787 IF( orgati ) THEN
788 a = z( i )*z( i ) + delta( ip1 )*
789 $ delta( ip1 )*( dpsi+dphi )
790 ELSE
791 a = z( ip1 )*z( ip1 ) +
792 $ delta( i )*delta( i )*( dpsi+dphi )
793 END IF
794 ELSE
795 a = delta( i )*delta( i )*dpsi +
796 $ delta( ip1 )*delta( ip1 )*dphi
797 END IF
798 END IF
799 eta = b / a
800 ELSE IF( a.LE.zero ) THEN
801 eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
802 ELSE
803 eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
804 END IF
805 ELSE
806*
807* Interpolation using THREE most relevant poles
808*
809 temp = rhoinv + psi + phi
810 IF( swtch ) THEN
811 c = temp - delta( iim1 )*dpsi - delta( iip1 )*dphi
812 zz( 1 ) = delta( iim1 )*delta( iim1 )*dpsi
813 zz( 3 ) = delta( iip1 )*delta( iip1 )*dphi
814 ELSE
815 IF( orgati ) THEN
816 temp1 = z( iim1 ) / delta( iim1 )
817 temp1 = temp1*temp1
818 c = temp - delta( iip1 )*( dpsi+dphi ) -
819 $ ( d( iim1 )-d( iip1 ) )*temp1
820 zz( 1 ) = z( iim1 )*z( iim1 )
821 zz( 3 ) = delta( iip1 )*delta( iip1 )*
822 $ ( ( dpsi-temp1 )+dphi )
823 ELSE
824 temp1 = z( iip1 ) / delta( iip1 )
825 temp1 = temp1*temp1
826 c = temp - delta( iim1 )*( dpsi+dphi ) -
827 $ ( d( iip1 )-d( iim1 ) )*temp1
828 zz( 1 ) = delta( iim1 )*delta( iim1 )*
829 $ ( dpsi+( dphi-temp1 ) )
830 zz( 3 ) = z( iip1 )*z( iip1 )
831 END IF
832 END IF
833 CALL dlaed6( niter, orgati, c, delta( iim1 ), zz, w,
834 $ eta,
835 $ info )
836 IF( info.NE.0 )
837 $ GO TO 250
838 END IF
839*
840* Note, eta should be positive if w is negative, and
841* eta should be negative otherwise. However,
842* if for some reason caused by roundoff, eta*w > 0,
843* we simply use one Newton step instead. This way
844* will guarantee eta*w < 0.
845*
846 IF( w*eta.GE.zero )
847 $ eta = -w / dw
848 temp = tau + eta
849 IF( temp.GT.dltub .OR. temp.LT.dltlb ) THEN
850 IF( w.LT.zero ) THEN
851 eta = ( dltub-tau ) / two
852 ELSE
853 eta = ( dltlb-tau ) / two
854 END IF
855 END IF
856*
857 DO 210 j = 1, n
858 delta( j ) = delta( j ) - eta
859 210 CONTINUE
860*
861 tau = tau + eta
862 prew = w
863*
864* Evaluate PSI and the derivative DPSI
865*
866 dpsi = zero
867 psi = zero
868 erretm = zero
869 DO 220 j = 1, iim1
870 temp = z( j ) / delta( j )
871 psi = psi + z( j )*temp
872 dpsi = dpsi + temp*temp
873 erretm = erretm + psi
874 220 CONTINUE
875 erretm = abs( erretm )
876*
877* Evaluate PHI and the derivative DPHI
878*
879 dphi = zero
880 phi = zero
881 DO 230 j = n, iip1, -1
882 temp = z( j ) / delta( j )
883 phi = phi + z( j )*temp
884 dphi = dphi + temp*temp
885 erretm = erretm + phi
886 230 CONTINUE
887*
888 temp = z( ii ) / delta( ii )
889 dw = dpsi + dphi + temp*temp
890 temp = z( ii )*temp
891 w = rhoinv + phi + psi + temp
892 erretm = eight*( phi-psi ) + erretm + two*rhoinv +
893 $ three*abs( temp ) + abs( tau )*dw
894 IF( w*prew.GT.zero .AND. abs( w ).GT.abs( prew ) / ten )
895 $ swtch = .NOT.swtch
896*
897 240 CONTINUE
898*
899* Return with INFO = 1, NITER = MAXIT and not converged
900*
901 info = 1
902 IF( orgati ) THEN
903 dlam = d( i ) + tau
904 ELSE
905 dlam = d( ip1 ) + tau
906 END IF
907*
908 END IF
909*
910 250 CONTINUE
911*
912 RETURN
913*
914* End of DLAED4
915*
916 END
dlaed4
subroutine dlaed4(n, i, d, z, delta, rho, dlam, info)
DLAED4 used by DSTEDC. Finds a single root of the secular equation.
Definition dlaed4.f:143
dlaed5
subroutine dlaed5(i, d, z, delta, rho, dlam)
DLAED5 used by DSTEDC. Solves the 2-by-2 secular equation.
Definition dlaed5.f:106
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:139