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 sbdsdc.
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 = 64 )
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
181  INTEGER ii, iim1, iip1, ip1, iter, j, niter
182  DOUBLE PRECISION a, b, c, delsq, delsq2, dphi, dpsi, dtiim,
183  $ dtiip, dtipsq, dtisq, dtnsq, dtnsq1, dw, eps,
184  $ erretm, eta, phi, prew, psi, rhoinv, sg2lb,
185  $ sg2ub, tau, 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 TAU 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  tau = two*b / ( sqrt( a*a+four*b*c )-a )
275  ELSE
276  tau = ( 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+TAU <= 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 TAU is to approximate
289 * SIGMA_n^2 - D( N )*D( N )
290 *
291  IF( a.LT.zero ) THEN
292  tau = two*b / ( sqrt( a*a+four*b*c )-a )
293  ELSE
294  tau = ( 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+TAU < SIGMA(N)^2 < D(N)^2+RHO/2
299 *
300  END IF
301 *
302 * The following ETA is to approximate SIGMA_n - D( N )
303 *
304  eta = tau / ( d( n )+sqrt( d( n )*d( n )+tau ) )
305 *
306  sigma = d( n ) + eta
307  DO 30 j = 1, n
308  delta( j ) = ( d( j )-d( i ) ) - eta
309  work( j ) = d( j ) + d( i ) + eta
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( tau )*( 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  tau = tau + eta
372  eta = eta / ( sigma+sqrt( eta+sigma*sigma ) )
373  DO 50 j = 1, n
374  delta( j ) = delta( j ) - eta
375  work( j ) = work( j ) + eta
376  50 continue
377 *
378  sigma = sigma + eta
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  temp = z( n ) / ( work( n )*delta( n ) )
396  phi = z( n )*temp
397  dphi = temp*temp
398  erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
399  $ abs( tau )*( dpsi+dphi )
400 *
401  w = rhoinv + phi + psi
402 *
403 * Main loop to update the values of the array DELTA
404 *
405  iter = niter + 1
406 *
407  DO 90 niter = iter, maxit
408 *
409 * Test for convergence
410 *
411  IF( abs( w ).LE.eps*erretm ) THEN
412  go to 240
413  END IF
414 *
415 * Calculate the new step
416 *
417  dtnsq1 = work( n-1 )*delta( n-1 )
418  dtnsq = work( n )*delta( n )
419  c = w - dtnsq1*dpsi - dtnsq*dphi
420  a = ( dtnsq+dtnsq1 )*w - dtnsq1*dtnsq*( dpsi+dphi )
421  b = dtnsq1*dtnsq*w
422  IF( a.GE.zero ) THEN
423  eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
424  ELSE
425  eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
426  END IF
427 *
428 * Note, eta should be positive if w is negative, and
429 * eta should be negative otherwise. However,
430 * if for some reason caused by roundoff, eta*w > 0,
431 * we simply use one Newton step instead. This way
432 * will guarantee eta*w < 0.
433 *
434  IF( w*eta.GT.zero )
435  $ eta = -w / ( dpsi+dphi )
436  temp = eta - dtnsq
437  IF( temp.LE.zero )
438  $ eta = eta / two
439 *
440  tau = tau + eta
441  eta = eta / ( sigma+sqrt( eta+sigma*sigma ) )
442  DO 70 j = 1, n
443  delta( j ) = delta( j ) - eta
444  work( j ) = work( j ) + eta
445  70 continue
446 *
447  sigma = sigma + eta
448 *
449 * Evaluate PSI and the derivative DPSI
450 *
451  dpsi = zero
452  psi = zero
453  erretm = zero
454  DO 80 j = 1, ii
455  temp = z( j ) / ( work( j )*delta( j ) )
456  psi = psi + z( j )*temp
457  dpsi = dpsi + temp*temp
458  erretm = erretm + psi
459  80 continue
460  erretm = abs( erretm )
461 *
462 * Evaluate PHI and the derivative DPHI
463 *
464  temp = z( n ) / ( work( n )*delta( n ) )
465  phi = z( n )*temp
466  dphi = temp*temp
467  erretm = eight*( -phi-psi ) + erretm - phi + rhoinv +
468  $ abs( tau )*( dpsi+dphi )
469 *
470  w = rhoinv + phi + psi
471  90 continue
472 *
473 * Return with INFO = 1, NITER = MAXIT and not converged
474 *
475  info = 1
476  go to 240
477 *
478 * End for the case I = N
479 *
480  ELSE
481 *
482 * The case for I < N
483 *
484  niter = 1
485  ip1 = i + 1
486 *
487 * Calculate initial guess
488 *
489  delsq = ( d( ip1 )-d( i ) )*( d( ip1 )+d( i ) )
490  delsq2 = delsq / two
491  temp = delsq2 / ( d( i )+sqrt( d( i )*d( i )+delsq2 ) )
492  DO 100 j = 1, n
493  work( j ) = d( j ) + d( i ) + temp
494  delta( j ) = ( d( j )-d( i ) ) - temp
495  100 continue
496 *
497  psi = zero
498  DO 110 j = 1, i - 1
499  psi = psi + z( j )*z( j ) / ( work( j )*delta( j ) )
500  110 continue
501 *
502  phi = zero
503  DO 120 j = n, i + 2, -1
504  phi = phi + z( j )*z( j ) / ( work( j )*delta( j ) )
505  120 continue
506  c = rhoinv + psi + phi
507  w = c + z( i )*z( i ) / ( work( i )*delta( i ) ) +
508  $ z( ip1 )*z( ip1 ) / ( work( ip1 )*delta( ip1 ) )
509 *
510  IF( w.GT.zero ) THEN
511 *
512 * d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2
513 *
514 * We choose d(i) as origin.
515 *
516  orgati = .true.
517  sg2lb = zero
518  sg2ub = delsq2
519  a = c*delsq + z( i )*z( i ) + z( ip1 )*z( ip1 )
520  b = z( i )*z( i )*delsq
521  IF( a.GT.zero ) THEN
522  tau = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
523  ELSE
524  tau = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
525  END IF
526 *
527 * TAU now is an estimation of SIGMA^2 - D( I )^2. The
528 * following, however, is the corresponding estimation of
529 * SIGMA - D( I ).
530 *
531  eta = tau / ( d( i )+sqrt( d( i )*d( i )+tau ) )
532  ELSE
533 *
534 * (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2
535 *
536 * We choose d(i+1) as origin.
537 *
538  orgati = .false.
539  sg2lb = -delsq2
540  sg2ub = zero
541  a = c*delsq - z( i )*z( i ) - z( ip1 )*z( ip1 )
542  b = z( ip1 )*z( ip1 )*delsq
543  IF( a.LT.zero ) THEN
544  tau = two*b / ( a-sqrt( abs( a*a+four*b*c ) ) )
545  ELSE
546  tau = -( a+sqrt( abs( a*a+four*b*c ) ) ) / ( two*c )
547  END IF
548 *
549 * TAU now is an estimation of SIGMA^2 - D( IP1 )^2. The
550 * following, however, is the corresponding estimation of
551 * SIGMA - D( IP1 ).
552 *
553  eta = tau / ( d( ip1 )+sqrt( abs( d( ip1 )*d( ip1 )+
554  $ tau ) ) )
555  END IF
556 *
557  IF( orgati ) THEN
558  ii = i
559  sigma = d( i ) + eta
560  DO 130 j = 1, n
561  work( j ) = d( j ) + d( i ) + eta
562  delta( j ) = ( d( j )-d( i ) ) - eta
563  130 continue
564  ELSE
565  ii = i + 1
566  sigma = d( ip1 ) + eta
567  DO 140 j = 1, n
568  work( j ) = d( j ) + d( ip1 ) + eta
569  delta( j ) = ( d( j )-d( ip1 ) ) - eta
570  140 continue
571  END IF
572  iim1 = ii - 1
573  iip1 = ii + 1
574 *
575 * Evaluate PSI and the derivative DPSI
576 *
577  dpsi = zero
578  psi = zero
579  erretm = zero
580  DO 150 j = 1, iim1
581  temp = z( j ) / ( work( j )*delta( j ) )
582  psi = psi + z( j )*temp
583  dpsi = dpsi + temp*temp
584  erretm = erretm + psi
585  150 continue
586  erretm = abs( erretm )
587 *
588 * Evaluate PHI and the derivative DPHI
589 *
590  dphi = zero
591  phi = zero
592  DO 160 j = n, iip1, -1
593  temp = z( j ) / ( work( j )*delta( j ) )
594  phi = phi + z( j )*temp
595  dphi = dphi + temp*temp
596  erretm = erretm + phi
597  160 continue
598 *
599  w = rhoinv + phi + psi
600 *
601 * W is the value of the secular function with
602 * its ii-th element removed.
603 *
604  swtch3 = .false.
605  IF( orgati ) THEN
606  IF( w.LT.zero )
607  $ swtch3 = .true.
608  ELSE
609  IF( w.GT.zero )
610  $ swtch3 = .true.
611  END IF
612  IF( ii.EQ.1 .OR. ii.EQ.n )
613  $ swtch3 = .false.
614 *
615  temp = z( ii ) / ( work( ii )*delta( ii ) )
616  dw = dpsi + dphi + temp*temp
617  temp = z( ii )*temp
618  w = w + temp
619  erretm = eight*( phi-psi ) + erretm + two*rhoinv +
620  $ three*abs( temp ) + abs( tau )*dw
621 *
622 * Test for convergence
623 *
624  IF( abs( w ).LE.eps*erretm ) THEN
625  go to 240
626  END IF
627 *
628  IF( w.LE.zero ) THEN
629  sg2lb = max( sg2lb, tau )
630  ELSE
631  sg2ub = min( sg2ub, tau )
632  END IF
633 *
634 * Calculate the new step
635 *
636  niter = niter + 1
637  IF( .NOT.swtch3 ) THEN
638  dtipsq = work( ip1 )*delta( ip1 )
639  dtisq = work( i )*delta( i )
640  IF( orgati ) THEN
641  c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
642  ELSE
643  c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
644  END IF
645  a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
646  b = dtipsq*dtisq*w
647  IF( c.EQ.zero ) THEN
648  IF( a.EQ.zero ) THEN
649  IF( orgati ) THEN
650  a = z( i )*z( i ) + dtipsq*dtipsq*( dpsi+dphi )
651  ELSE
652  a = z( ip1 )*z( ip1 ) + dtisq*dtisq*( dpsi+dphi )
653  END IF
654  END IF
655  eta = b / a
656  ELSE IF( a.LE.zero ) THEN
657  eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
658  ELSE
659  eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
660  END IF
661  ELSE
662 *
663 * Interpolation using THREE most relevant poles
664 *
665  dtiim = work( iim1 )*delta( iim1 )
666  dtiip = work( iip1 )*delta( iip1 )
667  temp = rhoinv + psi + phi
668  IF( orgati ) THEN
669  temp1 = z( iim1 ) / dtiim
670  temp1 = temp1*temp1
671  c = ( temp - dtiip*( dpsi+dphi ) ) -
672  $ ( d( iim1 )-d( iip1 ) )*( d( iim1 )+d( iip1 ) )*temp1
673  zz( 1 ) = z( iim1 )*z( iim1 )
674  IF( dpsi.LT.temp1 ) THEN
675  zz( 3 ) = dtiip*dtiip*dphi
676  ELSE
677  zz( 3 ) = dtiip*dtiip*( ( dpsi-temp1 )+dphi )
678  END IF
679  ELSE
680  temp1 = z( iip1 ) / dtiip
681  temp1 = temp1*temp1
682  c = ( temp - dtiim*( dpsi+dphi ) ) -
683  $ ( d( iip1 )-d( iim1 ) )*( d( iim1 )+d( iip1 ) )*temp1
684  IF( dphi.LT.temp1 ) THEN
685  zz( 1 ) = dtiim*dtiim*dpsi
686  ELSE
687  zz( 1 ) = dtiim*dtiim*( dpsi+( dphi-temp1 ) )
688  END IF
689  zz( 3 ) = z( iip1 )*z( iip1 )
690  END IF
691  zz( 2 ) = z( ii )*z( ii )
692  dd( 1 ) = dtiim
693  dd( 2 ) = delta( ii )*work( ii )
694  dd( 3 ) = dtiip
695  CALL dlaed6( niter, orgati, c, dd, zz, w, eta, info )
696  IF( info.NE.0 )
697  $ go to 240
698  END IF
699 *
700 * Note, eta should be positive if w is negative, and
701 * eta should be negative otherwise. However,
702 * if for some reason caused by roundoff, eta*w > 0,
703 * we simply use one Newton step instead. This way
704 * will guarantee eta*w < 0.
705 *
706  IF( w*eta.GE.zero )
707  $ eta = -w / dw
708  IF( orgati ) THEN
709  temp1 = work( i )*delta( i )
710  temp = eta - temp1
711  ELSE
712  temp1 = work( ip1 )*delta( ip1 )
713  temp = eta - temp1
714  END IF
715  IF( temp.GT.sg2ub .OR. temp.LT.sg2lb ) THEN
716  IF( w.LT.zero ) THEN
717  eta = ( sg2ub-tau ) / two
718  ELSE
719  eta = ( sg2lb-tau ) / two
720  END IF
721  END IF
722 *
723  tau = tau + eta
724  eta = eta / ( sigma+sqrt( sigma*sigma+eta ) )
725 *
726  prew = w
727 *
728  sigma = sigma + eta
729  DO 170 j = 1, n
730  work( j ) = work( j ) + eta
731  delta( j ) = delta( j ) - eta
732  170 continue
733 *
734 * Evaluate PSI and the derivative DPSI
735 *
736  dpsi = zero
737  psi = zero
738  erretm = zero
739  DO 180 j = 1, iim1
740  temp = z( j ) / ( work( j )*delta( j ) )
741  psi = psi + z( j )*temp
742  dpsi = dpsi + temp*temp
743  erretm = erretm + psi
744  180 continue
745  erretm = abs( erretm )
746 *
747 * Evaluate PHI and the derivative DPHI
748 *
749  dphi = zero
750  phi = zero
751  DO 190 j = n, iip1, -1
752  temp = z( j ) / ( work( j )*delta( j ) )
753  phi = phi + z( j )*temp
754  dphi = dphi + temp*temp
755  erretm = erretm + phi
756  190 continue
757 *
758  temp = z( ii ) / ( work( ii )*delta( ii ) )
759  dw = dpsi + dphi + temp*temp
760  temp = z( ii )*temp
761  w = rhoinv + phi + psi + temp
762  erretm = eight*( phi-psi ) + erretm + two*rhoinv +
763  $ three*abs( temp ) + abs( tau )*dw
764 *
765  IF( w.LE.zero ) THEN
766  sg2lb = max( sg2lb, tau )
767  ELSE
768  sg2ub = min( sg2ub, tau )
769  END IF
770 *
771  swtch = .false.
772  IF( orgati ) THEN
773  IF( -w.GT.abs( prew ) / ten )
774  $ swtch = .true.
775  ELSE
776  IF( w.GT.abs( prew ) / ten )
777  $ swtch = .true.
778  END IF
779 *
780 * Main loop to update the values of the array DELTA and WORK
781 *
782  iter = niter + 1
783 *
784  DO 230 niter = iter, maxit
785 *
786 * Test for convergence
787 *
788  IF( abs( w ).LE.eps*erretm ) THEN
789  go to 240
790  END IF
791 *
792 * Calculate the new step
793 *
794  IF( .NOT.swtch3 ) THEN
795  dtipsq = work( ip1 )*delta( ip1 )
796  dtisq = work( i )*delta( i )
797  IF( .NOT.swtch ) THEN
798  IF( orgati ) THEN
799  c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
800  ELSE
801  c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
802  END IF
803  ELSE
804  temp = z( ii ) / ( work( ii )*delta( ii ) )
805  IF( orgati ) THEN
806  dpsi = dpsi + temp*temp
807  ELSE
808  dphi = dphi + temp*temp
809  END IF
810  c = w - dtisq*dpsi - dtipsq*dphi
811  END IF
812  a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
813  b = dtipsq*dtisq*w
814  IF( c.EQ.zero ) THEN
815  IF( a.EQ.zero ) THEN
816  IF( .NOT.swtch ) THEN
817  IF( orgati ) THEN
818  a = z( i )*z( i ) + dtipsq*dtipsq*
819  $ ( dpsi+dphi )
820  ELSE
821  a = z( ip1 )*z( ip1 ) +
822  $ dtisq*dtisq*( dpsi+dphi )
823  END IF
824  ELSE
825  a = dtisq*dtisq*dpsi + dtipsq*dtipsq*dphi
826  END IF
827  END IF
828  eta = b / a
829  ELSE IF( a.LE.zero ) THEN
830  eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
831  ELSE
832  eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
833  END IF
834  ELSE
835 *
836 * Interpolation using THREE most relevant poles
837 *
838  dtiim = work( iim1 )*delta( iim1 )
839  dtiip = work( iip1 )*delta( iip1 )
840  temp = rhoinv + psi + phi
841  IF( swtch ) THEN
842  c = temp - dtiim*dpsi - dtiip*dphi
843  zz( 1 ) = dtiim*dtiim*dpsi
844  zz( 3 ) = dtiip*dtiip*dphi
845  ELSE
846  IF( orgati ) THEN
847  temp1 = z( iim1 ) / dtiim
848  temp1 = temp1*temp1
849  temp2 = ( d( iim1 )-d( iip1 ) )*
850  $ ( d( iim1 )+d( iip1 ) )*temp1
851  c = temp - dtiip*( dpsi+dphi ) - temp2
852  zz( 1 ) = z( iim1 )*z( iim1 )
853  IF( dpsi.LT.temp1 ) THEN
854  zz( 3 ) = dtiip*dtiip*dphi
855  ELSE
856  zz( 3 ) = dtiip*dtiip*( ( dpsi-temp1 )+dphi )
857  END IF
858  ELSE
859  temp1 = z( iip1 ) / dtiip
860  temp1 = temp1*temp1
861  temp2 = ( d( iip1 )-d( iim1 ) )*
862  $ ( d( iim1 )+d( iip1 ) )*temp1
863  c = temp - dtiim*( dpsi+dphi ) - temp2
864  IF( dphi.LT.temp1 ) THEN
865  zz( 1 ) = dtiim*dtiim*dpsi
866  ELSE
867  zz( 1 ) = dtiim*dtiim*( dpsi+( dphi-temp1 ) )
868  END IF
869  zz( 3 ) = z( iip1 )*z( iip1 )
870  END IF
871  END IF
872  dd( 1 ) = dtiim
873  dd( 2 ) = delta( ii )*work( ii )
874  dd( 3 ) = dtiip
875  CALL dlaed6( niter, orgati, c, dd, zz, w, eta, info )
876  IF( info.NE.0 )
877  $ go to 240
878  END IF
879 *
880 * Note, eta should be positive if w is negative, and
881 * eta should be negative otherwise. However,
882 * if for some reason caused by roundoff, eta*w > 0,
883 * we simply use one Newton step instead. This way
884 * will guarantee eta*w < 0.
885 *
886  IF( w*eta.GE.zero )
887  $ eta = -w / dw
888  IF( orgati ) THEN
889  temp1 = work( i )*delta( i )
890  temp = eta - temp1
891  ELSE
892  temp1 = work( ip1 )*delta( ip1 )
893  temp = eta - temp1
894  END IF
895  IF( temp.GT.sg2ub .OR. temp.LT.sg2lb ) THEN
896  IF( w.LT.zero ) THEN
897  eta = ( sg2ub-tau ) / two
898  ELSE
899  eta = ( sg2lb-tau ) / two
900  END IF
901  END IF
902 *
903  tau = tau + eta
904  eta = eta / ( sigma+sqrt( sigma*sigma+eta ) )
905 *
906  sigma = sigma + eta
907  DO 200 j = 1, n
908  work( j ) = work( j ) + eta
909  delta( j ) = delta( j ) - eta
910  200 continue
911 *
912  prew = w
913 *
914 * Evaluate PSI and the derivative DPSI
915 *
916  dpsi = zero
917  psi = zero
918  erretm = zero
919  DO 210 j = 1, iim1
920  temp = z( j ) / ( work( j )*delta( j ) )
921  psi = psi + z( j )*temp
922  dpsi = dpsi + temp*temp
923  erretm = erretm + psi
924  210 continue
925  erretm = abs( erretm )
926 *
927 * Evaluate PHI and the derivative DPHI
928 *
929  dphi = zero
930  phi = zero
931  DO 220 j = n, iip1, -1
932  temp = z( j ) / ( work( j )*delta( j ) )
933  phi = phi + z( j )*temp
934  dphi = dphi + temp*temp
935  erretm = erretm + phi
936  220 continue
937 *
938  temp = z( ii ) / ( work( ii )*delta( ii ) )
939  dw = dpsi + dphi + temp*temp
940  temp = z( ii )*temp
941  w = rhoinv + phi + psi + temp
942  erretm = eight*( phi-psi ) + erretm + two*rhoinv +
943  $ three*abs( temp ) + abs( tau )*dw
944  IF( w*prew.GT.zero .AND. abs( w ).GT.abs( prew ) / ten )
945  $ swtch = .NOT.swtch
946 *
947  IF( w.LE.zero ) THEN
948  sg2lb = max( sg2lb, tau )
949  ELSE
950  sg2ub = min( sg2ub, tau )
951  END IF
952 *
953  230 continue
954 *
955 * Return with INFO = 1, NITER = MAXIT and not converged
956 *
957  info = 1
958 *
959  END IF
960 *
961  240 continue
962  return
963 *
964 * End of DLASD4
965 *
966  END