OpenRadioss 2025.1.11
OpenRadioss project
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 auxOTHERcomputational
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 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 dlaed6( 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 dlaed6( 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 DLAED4
913*
914 END
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
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
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21