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