LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlaed4.f
Go to the documentation of this file.
1 *> \brief \b DLAED4 used by sstedc. 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 .GT. 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 *> \date September 2012
136 *
137 *> \ingroup auxOTHERcomputational
138 *
139 *> \par Contributors:
140 * ==================
141 *>
142 *> Ren-Cang Li, Computer Science Division, University of California
143 *> at Berkeley, USA
144 *>
145 * =====================================================================
146  SUBROUTINE dlaed4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
147 *
148 * -- LAPACK computational routine (version 3.4.2) --
149 * -- LAPACK is a software package provided by Univ. of Tennessee, --
150 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151 * September 2012
152 *
153 * .. Scalar Arguments ..
154  INTEGER i, info, n
155  DOUBLE PRECISION dlam, rho
156 * ..
157 * .. Array Arguments ..
158  DOUBLE PRECISION d( * ), delta( * ), z( * )
159 * ..
160 *
161 * =====================================================================
162 *
163 * .. Parameters ..
164  INTEGER maxit
165  parameter( maxit = 30 )
166  DOUBLE PRECISION zero, one, two, three, four, eight, ten
167  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
168  $ three = 3.0d0, four = 4.0d0, eight = 8.0d0,
169  $ ten = 10.0d0 )
170 * ..
171 * .. Local Scalars ..
172  LOGICAL orgati, swtch, swtch3
173  INTEGER ii, iim1, iip1, ip1, iter, j, niter
174  DOUBLE PRECISION a, b, c, del, dltlb, dltub, dphi, dpsi, dw,
175  $ eps, erretm, eta, midpt, phi, prew, psi,
176  $ rhoinv, tau, temp, temp1, w
177 * ..
178 * .. Local Arrays ..
179  DOUBLE PRECISION zz( 3 )
180 * ..
181 * .. External Functions ..
182  DOUBLE PRECISION dlamch
183  EXTERNAL dlamch
184 * ..
185 * .. External Subroutines ..
186  EXTERNAL dlaed5, dlaed6
187 * ..
188 * .. Intrinsic Functions ..
189  INTRINSIC abs, max, min, sqrt
190 * ..
191 * .. Executable Statements ..
192 *
193 * Since this routine is called in an inner loop, we do no argument
194 * checking.
195 *
196 * Quick return for N=1 and 2.
197 *
198  info = 0
199  IF( n.EQ.1 ) THEN
200 *
201 * Presumably, I=1 upon entry
202 *
203  dlam = d( 1 ) + rho*z( 1 )*z( 1 )
204  delta( 1 ) = one
205  return
206  END IF
207  IF( n.EQ.2 ) THEN
208  CALL dlaed5( i, d, z, delta, rho, dlam )
209  return
210  END IF
211 *
212 * Compute machine epsilon
213 *
214  eps = dlamch( 'Epsilon' )
215  rhoinv = one / rho
216 *
217 * The case I = N
218 *
219  IF( i.EQ.n ) THEN
220 *
221 * Initialize some basic variables
222 *
223  ii = n - 1
224  niter = 1
225 *
226 * Calculate initial guess
227 *
228  midpt = rho / two
229 *
230 * If ||Z||_2 is not one, then TEMP should be set to
231 * RHO * ||Z||_2^2 / TWO
232 *
233  DO 10 j = 1, n
234  delta( j ) = ( d( j )-d( i ) ) - midpt
235  10 continue
236 *
237  psi = zero
238  DO 20 j = 1, n - 2
239  psi = psi + z( j )*z( j ) / delta( j )
240  20 continue
241 *
242  c = rhoinv + psi
243  w = c + z( ii )*z( ii ) / delta( ii ) +
244  $ z( n )*z( n ) / delta( n )
245 *
246  IF( w.LE.zero ) THEN
247  temp = z( n-1 )*z( n-1 ) / ( d( n )-d( n-1 )+rho ) +
248  $ z( n )*z( n ) / rho
249  IF( c.LE.temp ) THEN
250  tau = rho
251  ELSE
252  del = d( n ) - d( n-1 )
253  a = -c*del + z( n-1 )*z( n-1 ) + z( n )*z( n )
254  b = z( n )*z( n )*del
255  IF( a.LT.zero ) THEN
256  tau = two*b / ( sqrt( a*a+four*b*c )-a )
257  ELSE
258  tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
259  END IF
260  END IF
261 *
262 * It can be proved that
263 * D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO
264 *
265  dltlb = midpt
266  dltub = rho
267  ELSE
268  del = d( n ) - d( n-1 )
269  a = -c*del + z( n-1 )*z( n-1 ) + z( n )*z( n )
270  b = z( n )*z( n )*del
271  IF( a.LT.zero ) THEN
272  tau = two*b / ( sqrt( a*a+four*b*c )-a )
273  ELSE
274  tau = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
275  END IF
276 *
277 * It can be proved that
278 * D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2
279 *
280  dltlb = zero
281  dltub = midpt
282  END IF
283 *
284  DO 30 j = 1, n
285  delta( j ) = ( d( j )-d( i ) ) - tau
286  30 continue
287 *
288 * Evaluate PSI and the derivative DPSI
289 *
290  dpsi = zero
291  psi = zero
292  erretm = zero
293  DO 40 j = 1, ii
294  temp = z( j ) / delta( j )
295  psi = psi + z( j )*temp
296  dpsi = dpsi + temp*temp
297  erretm = erretm + psi
298  40 continue
299  erretm = abs( erretm )
300 *
301 * Evaluate PHI and the derivative DPHI
302 *
303  temp = z( n ) / delta( n )
304  phi = z( n )*temp
305  dphi = temp*temp
306  erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
307  $ abs( tau )*( dpsi+dphi )
308 *
309  w = rhoinv + phi + psi
310 *
311 * Test for convergence
312 *
313  IF( abs( w ).LE.eps*erretm ) THEN
314  dlam = d( i ) + tau
315  go to 250
316  END IF
317 *
318  IF( w.LE.zero ) THEN
319  dltlb = max( dltlb, tau )
320  ELSE
321  dltub = min( dltub, tau )
322  END IF
323 *
324 * Calculate the new step
325 *
326  niter = niter + 1
327  c = w - delta( n-1 )*dpsi - delta( n )*dphi
328  a = ( delta( n-1 )+delta( n ) )*w -
329  $ delta( n-1 )*delta( n )*( dpsi+dphi )
330  b = delta( n-1 )*delta( n )*w
331  IF( c.LT.zero )
332  $ c = abs( c )
333  IF( c.EQ.zero ) THEN
334 * ETA = B/A
335 * ETA = RHO - TAU
336  eta = dltub - tau
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