LAPACK 3.12.0
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*> \htmlonly
9*> Download DLAED4 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed4.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed4.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed4.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER I, INFO, N
25* DOUBLE PRECISION DLAM, RHO
26* ..
27* .. Array Arguments ..
28* DOUBLE PRECISION D( * ), DELTA( * ), Z( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> This subroutine computes the I-th updated eigenvalue of a symmetric
38*> rank-one modification to a diagonal matrix whose elements are
39*> given in the array d, and that
40*>
41*> D(i) < D(j) for i < j
42*>
43*> and that RHO > 0. This is arranged by the calling routine, and is
44*> no loss in generality. The rank-one modified system is thus
45*>
46*> diag( D ) + RHO * Z * Z_transpose.
47*>
48*> where we assume the Euclidean norm of Z is 1.
49*>
50*> The method consists of approximating the rational functions in the
51*> secular equation by simpler interpolating rational functions.
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] N
58*> \verbatim
59*> N is INTEGER
60*> The length of all arrays.
61*> \endverbatim
62*>
63*> \param[in] I
64*> \verbatim
65*> I is INTEGER
66*> The index of the eigenvalue to be computed. 1 <= I <= N.
67*> \endverbatim
68*>
69*> \param[in] D
70*> \verbatim
71*> D is DOUBLE PRECISION array, dimension (N)
72*> The original eigenvalues. It is assumed that they are in
73*> order, D(I) < D(J) for I < J.
74*> \endverbatim
75*>
76*> \param[in] Z
77*> \verbatim
78*> Z is DOUBLE PRECISION array, dimension (N)
79*> The components of the updating vector.
80*> \endverbatim
81*>
82*> \param[out] DELTA
83*> \verbatim
84*> DELTA is DOUBLE PRECISION array, dimension (N)
85*> If N > 2, DELTA contains (D(j) - lambda_I) in its j-th
86*> component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5
87*> for detail. The vector DELTA contains the information necessary
88*> to construct the eigenvectors by DLAED3 and DLAED9.
89*> \endverbatim
90*>
91*> \param[in] RHO
92*> \verbatim
93*> RHO is DOUBLE PRECISION
94*> The scalar in the symmetric updating formula.
95*> \endverbatim
96*>
97*> \param[out] DLAM
98*> \verbatim
99*> DLAM is DOUBLE PRECISION
100*> The computed lambda_I, the I-th updated eigenvalue.
101*> \endverbatim
102*>
103*> \param[out] INFO
104*> \verbatim
105*> INFO is INTEGER
106*> = 0: successful exit
107*> > 0: if INFO = 1, the updating process failed.
108*> \endverbatim
109*
110*> \par Internal Parameters:
111* =========================
112*>
113*> \verbatim
114*> Logical variable ORGATI (origin-at-i?) is used for distinguishing
115*> whether D(i) or D(i+1) is treated as the origin.
116*>
117*> ORGATI = .true. origin at i
118*> ORGATI = .false. origin at i+1
119*>
120*> Logical variable SWTCH3 (switch-for-3-poles?) is for noting
121*> if we are working with THREE poles!
122*>
123*> MAXIT is the maximum number of iterations allowed for each
124*> eigenvalue.
125*> \endverbatim
126*
127* Authors:
128* ========
129*
130*> \author Univ. of Tennessee
131*> \author Univ. of California Berkeley
132*> \author Univ. of Colorado Denver
133*> \author NAG Ltd.
134*
135*> \ingroup laed4
136*
137*> \par Contributors:
138* ==================
139*>
140*> Ren-Cang Li, Computer Science Division, University of California
141*> at Berkeley, USA
142*>
143* =====================================================================
144 SUBROUTINE dlaed4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
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*
917 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:145
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