LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slasd4.f
Go to the documentation of this file.
1 *> \brief \b SLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. Used by sbdsdc.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SLASD4 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd4.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd4.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd4.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
22 *
23 * .. Scalar Arguments ..
24 * INTEGER I, INFO, N
25 * REAL RHO, SIGMA
26 * ..
27 * .. Array Arguments ..
28 * REAL 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 REAL 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 REAL array, dimension ( N )
80 *> The components of the updating vector.
81 *> \endverbatim
82 *>
83 *> \param[out] DELTA
84 *> \verbatim
85 *> DELTA is REAL 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 REAL
95 *> The scalar in the symmetric updating formula.
96 *> \endverbatim
97 *>
98 *> \param[out] SIGMA
99 *> \verbatim
100 *> SIGMA is REAL
101 *> The computed sigma_I, the I-th updated eigenvalue.
102 *> \endverbatim
103 *>
104 *> \param[out] WORK
105 *> \verbatim
106 *> WORK is REAL 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 slasd4( 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  REAL rho, sigma
164 * ..
165 * .. Array Arguments ..
166  REAL d( * ), delta( * ), work( * ), z( * )
167 * ..
168 *
169 * =====================================================================
170 *
171 * .. Parameters ..
172  INTEGER maxit
173  parameter( maxit = 400 )
174  REAL zero, one, two, three, four, eight, ten
175  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0,
176  $ three = 3.0e+0, four = 4.0e+0, eight = 8.0e+0,
177  $ ten = 10.0e+0 )
178 * ..
179 * .. Local Scalars ..
180  LOGICAL orgati, swtch, swtch3, geomavg
181  INTEGER ii, iim1, iip1, ip1, iter, j, niter
182  REAL 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  REAL dd( 3 ), zz( 3 )
189 * ..
190 * .. External Subroutines ..
191  EXTERNAL slaed6, slasd5
192 * ..
193 * .. External Functions ..
194  REAL slamch
195  EXTERNAL slamch
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 slasd5( i, d, z, delta, rho, sigma, work )
219  RETURN
220  END IF
221 *
222 * Compute machine epsilon
223 *
224  eps = slamch( '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 slaed6( niter, orgati, c, dd, zz, w, eta, info )
699 *
700  IF( info.NE.0 ) THEN
701 *
702 * If INFO is not 0, i.e., SLAED6 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 slaed6( niter, orgati, c, dd, zz, w, eta, info )
918 *
919  IF( info.NE.0 ) THEN
920 *
921 * If INFO is not 0, i.e., SLAED6 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 SLASD4
1056 *
1057  END