LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 November 2013
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.5.0) --
157 * -- LAPACK is a software package provided by Univ. of Tennessee, --
158 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159 * November 2013
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  tau2= zero
227 *
228 * The case I = N
229 *
230  IF( i.EQ.n ) THEN
231 *
232 * Initialize some basic variables
233 *
234  ii = n - 1
235  niter = 1
236 *
237 * Calculate initial guess
238 *
239  temp = rho / two
240 *
241 * If ||Z||_2 is not one, then TEMP should be set to
242 * RHO * ||Z||_2^2 / TWO
243 *
244  temp1 = temp / ( d( n )+sqrt( d( n )*d( n )+temp ) )
245  DO 10 j = 1, n
246  work( j ) = d( j ) + d( n ) + temp1
247  delta( j ) = ( d( j )-d( n ) ) - temp1
248  10 CONTINUE
249 *
250  psi = zero
251  DO 20 j = 1, n - 2
252  psi = psi + z( j )*z( j ) / ( delta( j )*work( j ) )
253  20 CONTINUE
254 *
255  c = rhoinv + psi
256  w = c + z( ii )*z( ii ) / ( delta( ii )*work( ii ) ) +
257  $ z( n )*z( n ) / ( delta( n )*work( n ) )
258 *
259  IF( w.LE.zero ) THEN
260  temp1 = sqrt( d( n )*d( n )+rho )
261  temp = z( n-1 )*z( n-1 ) / ( ( d( n-1 )+temp1 )*
262  $ ( d( n )-d( n-1 )+rho / ( d( n )+temp1 ) ) ) +
263  $ z( n )*z( n ) / rho
264 *
265 * The following TAU2 is to approximate
266 * SIGMA_n^2 - D( N )*D( N )
267 *
268  IF( c.LE.temp ) THEN
269  tau = rho
270  ELSE
271  delsq = ( d( n )-d( n-1 ) )*( d( n )+d( n-1 ) )
272  a = -c*delsq + z( n-1 )*z( n-1 ) + z( n )*z( n )
273  b = z( n )*z( n )*delsq
274  IF( a.LT.zero ) THEN
275  tau2 = two*b / ( sqrt( a*a+four*b*c )-a )
276  ELSE
277  tau2 = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
278  END IF
279  tau = tau2 / ( d( n )+sqrt( d( n )*d( n )+tau2 ) )
280  END IF
281 *
282 * It can be proved that
283 * D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU2 <= D(N)^2+RHO
284 *
285  ELSE
286  delsq = ( d( n )-d( n-1 ) )*( d( n )+d( n-1 ) )
287  a = -c*delsq + z( n-1 )*z( n-1 ) + z( n )*z( n )
288  b = z( n )*z( n )*delsq
289 *
290 * The following TAU2 is to approximate
291 * SIGMA_n^2 - D( N )*D( N )
292 *
293  IF( a.LT.zero ) THEN
294  tau2 = two*b / ( sqrt( a*a+four*b*c )-a )
295  ELSE
296  tau2 = ( a+sqrt( a*a+four*b*c ) ) / ( two*c )
297  END IF
298  tau = tau2 / ( d( n )+sqrt( d( n )*d( n )+tau2 ) )
299 
300 *
301 * It can be proved that
302 * D(N)^2 < D(N)^2+TAU2 < SIGMA(N)^2 < D(N)^2+RHO/2
303 *
304  END IF
305 *
306 * The following TAU is to approximate SIGMA_n - D( N )
307 *
308 * TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
309 *
310  sigma = d( n ) + tau
311  DO 30 j = 1, n
312  delta( j ) = ( d( j )-d( n ) ) - tau
313  work( j ) = d( j ) + d( n ) + tau
314  30 CONTINUE
315 *
316 * Evaluate PSI and the derivative DPSI
317 *
318  dpsi = zero
319  psi = zero
320  erretm = zero
321  DO 40 j = 1, ii
322  temp = z( j ) / ( delta( j )*work( j ) )
323  psi = psi + z( j )*temp
324  dpsi = dpsi + temp*temp
325  erretm = erretm + psi
326  40 CONTINUE
327  erretm = abs( erretm )
328 *
329 * Evaluate PHI and the derivative DPHI
330 *
331  temp = z( n ) / ( delta( n )*work( n ) )
332  phi = z( n )*temp
333  dphi = temp*temp
334  erretm = eight*( -phi-psi ) + erretm - phi + rhoinv
335 * $ + ABS( TAU2 )*( DPSI+DPHI )
336 *
337  w = rhoinv + phi + psi
338 *
339 * Test for convergence
340 *
341  IF( abs( w ).LE.eps*erretm ) THEN
342  GO TO 240
343  END IF
344 *
345 * Calculate the new step
346 *
347  niter = niter + 1
348  dtnsq1 = work( n-1 )*delta( n-1 )
349  dtnsq = work( n )*delta( n )
350  c = w - dtnsq1*dpsi - dtnsq*dphi
351  a = ( dtnsq+dtnsq1 )*w - dtnsq*dtnsq1*( dpsi+dphi )
352  b = dtnsq*dtnsq1*w
353  IF( c.LT.zero )
354  $ c = abs( c )
355  IF( c.EQ.zero ) THEN
356  eta = rho - sigma*sigma
357  ELSE IF( a.GE.zero ) THEN
358  eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
359  ELSE
360  eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
361  END IF
362 *
363 * Note, eta should be positive if w is negative, and
364 * eta should be negative otherwise. However,
365 * if for some reason caused by roundoff, eta*w > 0,
366 * we simply use one Newton step instead. This way
367 * will guarantee eta*w < 0.
368 *
369  IF( w*eta.GT.zero )
370  $ eta = -w / ( dpsi+dphi )
371  temp = eta - dtnsq
372  IF( temp.GT.rho )
373  $ eta = rho + dtnsq
374 *
375  eta = eta / ( sigma+sqrt( eta+sigma*sigma ) )
376  tau = tau + eta
377  sigma = sigma + eta
378 *
379  DO 50 j = 1, n
380  delta( j ) = delta( j ) - eta
381  work( j ) = work( j ) + eta
382  50 CONTINUE
383 *
384 * Evaluate PSI and the derivative DPSI
385 *
386  dpsi = zero
387  psi = zero
388  erretm = zero
389  DO 60 j = 1, ii
390  temp = z( j ) / ( work( j )*delta( j ) )
391  psi = psi + z( j )*temp
392  dpsi = dpsi + temp*temp
393  erretm = erretm + psi
394  60 CONTINUE
395  erretm = abs( erretm )
396 *
397 * Evaluate PHI and the derivative DPHI
398 *
399  tau2 = work( n )*delta( n )
400  temp = z( n ) / tau2
401  phi = z( n )*temp
402  dphi = temp*temp
403  erretm = eight*( -phi-psi ) + erretm - phi + rhoinv
404 * $ + ABS( TAU2 )*( DPSI+DPHI )
405 *
406  w = rhoinv + phi + psi
407 *
408 * Main loop to update the values of the array DELTA
409 *
410  iter = niter + 1
411 *
412  DO 90 niter = iter, maxit
413 *
414 * Test for convergence
415 *
416  IF( abs( w ).LE.eps*erretm ) THEN
417  GO TO 240
418  END IF
419 *
420 * Calculate the new step
421 *
422  dtnsq1 = work( n-1 )*delta( n-1 )
423  dtnsq = work( n )*delta( n )
424  c = w - dtnsq1*dpsi - dtnsq*dphi
425  a = ( dtnsq+dtnsq1 )*w - dtnsq1*dtnsq*( dpsi+dphi )
426  b = dtnsq1*dtnsq*w
427  IF( a.GE.zero ) THEN
428  eta = ( a+sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
429  ELSE
430  eta = two*b / ( a-sqrt( abs( a*a-four*b*c ) ) )
431  END IF
432 *
433 * Note, eta should be positive if w is negative, and
434 * eta should be negative otherwise. However,
435 * if for some reason caused by roundoff, eta*w > 0,
436 * we simply use one Newton step instead. This way
437 * will guarantee eta*w < 0.
438 *
439  IF( w*eta.GT.zero )
440  $ eta = -w / ( dpsi+dphi )
441  temp = eta - dtnsq
442  IF( temp.LE.zero )
443  $ eta = eta / two
444 *
445  eta = eta / ( sigma+sqrt( eta+sigma*sigma ) )
446  tau = tau + eta
447  sigma = sigma + eta
448 *
449  DO 70 j = 1, n
450  delta( j ) = delta( j ) - eta
451  work( j ) = work( j ) + eta
452  70 CONTINUE
453 *
454 * Evaluate PSI and the derivative DPSI
455 *
456  dpsi = zero
457  psi = zero
458  erretm = zero
459  DO 80 j = 1, ii
460  temp = z( j ) / ( work( j )*delta( j ) )
461  psi = psi + z( j )*temp
462  dpsi = dpsi + temp*temp
463  erretm = erretm + psi
464  80 CONTINUE
465  erretm = abs( erretm )
466 *
467 * Evaluate PHI and the derivative DPHI
468 *
469  tau2 = work( n )*delta( n )
470  temp = z( n ) / tau2
471  phi = z( n )*temp
472  dphi = temp*temp
473  erretm = eight*( -phi-psi ) + erretm - phi + rhoinv
474 * $ + ABS( TAU2 )*( DPSI+DPHI )
475 *
476  w = rhoinv + phi + psi
477  90 CONTINUE
478 *
479 * Return with INFO = 1, NITER = MAXIT and not converged
480 *
481  info = 1
482  GO TO 240
483 *
484 * End for the case I = N
485 *
486  ELSE
487 *
488 * The case for I < N
489 *
490  niter = 1
491  ip1 = i + 1
492 *
493 * Calculate initial guess
494 *
495  delsq = ( d( ip1 )-d( i ) )*( d( ip1 )+d( i ) )
496  delsq2 = delsq / two
497  sq2=sqrt( ( d( i )*d( i )+d( ip1 )*d( ip1 ) ) / two )
498  temp = delsq2 / ( d( i )+sq2 )
499  DO 100 j = 1, n
500  work( j ) = d( j ) + d( i ) + temp
501  delta( j ) = ( d( j )-d( i ) ) - temp
502  100 CONTINUE
503 *
504  psi = zero
505  DO 110 j = 1, i - 1
506  psi = psi + z( j )*z( j ) / ( work( j )*delta( j ) )
507  110 CONTINUE
508 *
509  phi = zero
510  DO 120 j = n, i + 2, -1
511  phi = phi + z( j )*z( j ) / ( work( j )*delta( j ) )
512  120 CONTINUE
513  c = rhoinv + psi + phi
514  w = c + z( i )*z( i ) / ( work( i )*delta( i ) ) +
515  $ z( ip1 )*z( ip1 ) / ( work( ip1 )*delta( ip1 ) )
516 *
517  geomavg = .false.
518  IF( w.GT.zero ) THEN
519 *
520 * d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2
521 *
522 * We choose d(i) as origin.
523 *
524  orgati = .true.
525  ii = i
526  sglb = zero
527  sgub = delsq2 / ( d( i )+sq2 )
528  a = c*delsq + z( i )*z( i ) + z( ip1 )*z( ip1 )
529  b = z( i )*z( i )*delsq
530  IF( a.GT.zero ) THEN
531  tau2 = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
532  ELSE
533  tau2 = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
534  END IF
535 *
536 * TAU2 now is an estimation of SIGMA^2 - D( I )^2. The
537 * following, however, is the corresponding estimation of
538 * SIGMA - D( I ).
539 *
540  tau = tau2 / ( d( i )+sqrt( d( i )*d( i )+tau2 ) )
541  temp = sqrt(eps)
542  IF( (d(i).LE.temp*d(ip1)).AND.(abs(z(i)).LE.temp)
543  $ .AND.(d(i).GT.zero) ) THEN
544  tau = min( ten*d(i), sgub )
545  geomavg = .true.
546  END IF
547  ELSE
548 *
549 * (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2
550 *
551 * We choose d(i+1) as origin.
552 *
553  orgati = .false.
554  ii = ip1
555  sglb = -delsq2 / ( d( ii )+sq2 )
556  sgub = zero
557  a = c*delsq - z( i )*z( i ) - z( ip1 )*z( ip1 )
558  b = z( ip1 )*z( ip1 )*delsq
559  IF( a.LT.zero ) THEN
560  tau2 = two*b / ( a-sqrt( abs( a*a+four*b*c ) ) )
561  ELSE
562  tau2 = -( a+sqrt( abs( a*a+four*b*c ) ) ) / ( two*c )
563  END IF
564 *
565 * TAU2 now is an estimation of SIGMA^2 - D( IP1 )^2. The
566 * following, however, is the corresponding estimation of
567 * SIGMA - D( IP1 ).
568 *
569  tau = tau2 / ( d( ip1 )+sqrt( abs( d( ip1 )*d( ip1 )+
570  $ tau2 ) ) )
571  END IF
572 *
573  sigma = d( ii ) + tau
574  DO 130 j = 1, n
575  work( j ) = d( j ) + d( ii ) + tau
576  delta( j ) = ( d( j )-d( ii ) ) - tau
577  130 CONTINUE
578  iim1 = ii - 1
579  iip1 = ii + 1
580 *
581 * Evaluate PSI and the derivative DPSI
582 *
583  dpsi = zero
584  psi = zero
585  erretm = zero
586  DO 150 j = 1, iim1
587  temp = z( j ) / ( work( j )*delta( j ) )
588  psi = psi + z( j )*temp
589  dpsi = dpsi + temp*temp
590  erretm = erretm + psi
591  150 CONTINUE
592  erretm = abs( erretm )
593 *
594 * Evaluate PHI and the derivative DPHI
595 *
596  dphi = zero
597  phi = zero
598  DO 160 j = n, iip1, -1
599  temp = z( j ) / ( work( j )*delta( j ) )
600  phi = phi + z( j )*temp
601  dphi = dphi + temp*temp
602  erretm = erretm + phi
603  160 CONTINUE
604 *
605  w = rhoinv + phi + psi
606 *
607 * W is the value of the secular function with
608 * its ii-th element removed.
609 *
610  swtch3 = .false.
611  IF( orgati ) THEN
612  IF( w.LT.zero )
613  $ swtch3 = .true.
614  ELSE
615  IF( w.GT.zero )
616  $ swtch3 = .true.
617  END IF
618  IF( ii.EQ.1 .OR. ii.EQ.n )
619  $ swtch3 = .false.
620 *
621  temp = z( ii ) / ( work( ii )*delta( ii ) )
622  dw = dpsi + dphi + temp*temp
623  temp = z( ii )*temp
624  w = w + temp
625  erretm = eight*( phi-psi ) + erretm + two*rhoinv
626  $ + three*abs( temp )
627 * $ + ABS( TAU2 )*DW
628 *
629 * Test for convergence
630 *
631  IF( abs( w ).LE.eps*erretm ) THEN
632  GO TO 240
633  END IF
634 *
635  IF( w.LE.zero ) THEN
636  sglb = max( sglb, tau )
637  ELSE
638  sgub = min( sgub, tau )
639  END IF
640 *
641 * Calculate the new step
642 *
643  niter = niter + 1
644  IF( .NOT.swtch3 ) THEN
645  dtipsq = work( ip1 )*delta( ip1 )
646  dtisq = work( i )*delta( i )
647  IF( orgati ) THEN
648  c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
649  ELSE
650  c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
651  END IF
652  a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
653  b = dtipsq*dtisq*w
654  IF( c.EQ.zero ) THEN
655  IF( a.EQ.zero ) THEN
656  IF( orgati ) THEN
657  a = z( i )*z( i ) + dtipsq*dtipsq*( dpsi+dphi )
658  ELSE
659  a = z( ip1 )*z( ip1 ) + dtisq*dtisq*( dpsi+dphi )
660  END IF
661  END IF
662  eta = b / a
663  ELSE IF( a.LE.zero ) THEN
664  eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
665  ELSE
666  eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
667  END IF
668  ELSE
669 *
670 * Interpolation using THREE most relevant poles
671 *
672  dtiim = work( iim1 )*delta( iim1 )
673  dtiip = work( iip1 )*delta( iip1 )
674  temp = rhoinv + psi + phi
675  IF( orgati ) THEN
676  temp1 = z( iim1 ) / dtiim
677  temp1 = temp1*temp1
678  c = ( temp - dtiip*( dpsi+dphi ) ) -
679  $ ( d( iim1 )-d( iip1 ) )*( d( iim1 )+d( iip1 ) )*temp1
680  zz( 1 ) = z( iim1 )*z( iim1 )
681  IF( dpsi.LT.temp1 ) THEN
682  zz( 3 ) = dtiip*dtiip*dphi
683  ELSE
684  zz( 3 ) = dtiip*dtiip*( ( dpsi-temp1 )+dphi )
685  END IF
686  ELSE
687  temp1 = z( iip1 ) / dtiip
688  temp1 = temp1*temp1
689  c = ( temp - dtiim*( dpsi+dphi ) ) -
690  $ ( d( iip1 )-d( iim1 ) )*( d( iim1 )+d( iip1 ) )*temp1
691  IF( dphi.LT.temp1 ) THEN
692  zz( 1 ) = dtiim*dtiim*dpsi
693  ELSE
694  zz( 1 ) = dtiim*dtiim*( dpsi+( dphi-temp1 ) )
695  END IF
696  zz( 3 ) = z( iip1 )*z( iip1 )
697  END IF
698  zz( 2 ) = z( ii )*z( ii )
699  dd( 1 ) = dtiim
700  dd( 2 ) = delta( ii )*work( ii )
701  dd( 3 ) = dtiip
702  CALL dlaed6( niter, orgati, c, dd, zz, w, eta, info )
703 *
704  IF( info.NE.0 ) THEN
705 *
706 * If INFO is not 0, i.e., DLAED6 failed, switch back
707 * to 2 pole interpolation.
708 *
709  swtch3 = .false.
710  info = 0
711  dtipsq = work( ip1 )*delta( ip1 )
712  dtisq = work( i )*delta( i )
713  IF( orgati ) THEN
714  c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
715  ELSE
716  c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
717  END IF
718  a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
719  b = dtipsq*dtisq*w
720  IF( c.EQ.zero ) THEN
721  IF( a.EQ.zero ) THEN
722  IF( orgati ) THEN
723  a = z( i )*z( i ) + dtipsq*dtipsq*( dpsi+dphi )
724  ELSE
725  a = z( ip1 )*z( ip1 ) + dtisq*dtisq*( dpsi+dphi)
726  END IF
727  END IF
728  eta = b / a
729  ELSE IF( a.LE.zero ) THEN
730  eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
731  ELSE
732  eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
733  END IF
734  END IF
735  END IF
736 *
737 * Note, eta should be positive if w is negative, and
738 * eta should be negative otherwise. However,
739 * if for some reason caused by roundoff, eta*w > 0,
740 * we simply use one Newton step instead. This way
741 * will guarantee eta*w < 0.
742 *
743  IF( w*eta.GE.zero )
744  $ eta = -w / dw
745 *
746  eta = eta / ( sigma+sqrt( sigma*sigma+eta ) )
747  temp = tau + eta
748  IF( temp.GT.sgub .OR. temp.LT.sglb ) THEN
749  IF( w.LT.zero ) THEN
750  eta = ( sgub-tau ) / two
751  ELSE
752  eta = ( sglb-tau ) / two
753  END IF
754  IF( geomavg ) THEN
755  IF( w .LT. zero ) THEN
756  IF( tau .GT. zero ) THEN
757  eta = sqrt(sgub*tau)-tau
758  END IF
759  ELSE
760  IF( sglb .GT. zero ) THEN
761  eta = sqrt(sglb*tau)-tau
762  END IF
763  END IF
764  END IF
765  END IF
766 *
767  prew = w
768 *
769  tau = tau + eta
770  sigma = sigma + eta
771 *
772  DO 170 j = 1, n
773  work( j ) = work( j ) + eta
774  delta( j ) = delta( j ) - eta
775  170 CONTINUE
776 *
777 * Evaluate PSI and the derivative DPSI
778 *
779  dpsi = zero
780  psi = zero
781  erretm = zero
782  DO 180 j = 1, iim1
783  temp = z( j ) / ( work( j )*delta( j ) )
784  psi = psi + z( j )*temp
785  dpsi = dpsi + temp*temp
786  erretm = erretm + psi
787  180 CONTINUE
788  erretm = abs( erretm )
789 *
790 * Evaluate PHI and the derivative DPHI
791 *
792  dphi = zero
793  phi = zero
794  DO 190 j = n, iip1, -1
795  temp = z( j ) / ( work( j )*delta( j ) )
796  phi = phi + z( j )*temp
797  dphi = dphi + temp*temp
798  erretm = erretm + phi
799  190 CONTINUE
800 *
801  tau2 = work( ii )*delta( ii )
802  temp = z( ii ) / tau2
803  dw = dpsi + dphi + temp*temp
804  temp = z( ii )*temp
805  w = rhoinv + phi + psi + temp
806  erretm = eight*( phi-psi ) + erretm + two*rhoinv
807  $ + three*abs( temp )
808 * $ + ABS( TAU2 )*DW
809 *
810  swtch = .false.
811  IF( orgati ) THEN
812  IF( -w.GT.abs( prew ) / ten )
813  $ swtch = .true.
814  ELSE
815  IF( w.GT.abs( prew ) / ten )
816  $ swtch = .true.
817  END IF
818 *
819 * Main loop to update the values of the array DELTA and WORK
820 *
821  iter = niter + 1
822 *
823  DO 230 niter = iter, maxit
824 *
825 * Test for convergence
826 *
827  IF( abs( w ).LE.eps*erretm ) THEN
828 * $ .OR. (SGUB-SGLB).LE.EIGHT*ABS(SGUB+SGLB) ) THEN
829  GO TO 240
830  END IF
831 *
832  IF( w.LE.zero ) THEN
833  sglb = max( sglb, tau )
834  ELSE
835  sgub = min( sgub, tau )
836  END IF
837 *
838 * Calculate the new step
839 *
840  IF( .NOT.swtch3 ) THEN
841  dtipsq = work( ip1 )*delta( ip1 )
842  dtisq = work( i )*delta( i )
843  IF( .NOT.swtch ) THEN
844  IF( orgati ) THEN
845  c = w - dtipsq*dw + delsq*( z( i ) / dtisq )**2
846  ELSE
847  c = w - dtisq*dw - delsq*( z( ip1 ) / dtipsq )**2
848  END IF
849  ELSE
850  temp = z( ii ) / ( work( ii )*delta( ii ) )
851  IF( orgati ) THEN
852  dpsi = dpsi + temp*temp
853  ELSE
854  dphi = dphi + temp*temp
855  END IF
856  c = w - dtisq*dpsi - dtipsq*dphi
857  END IF
858  a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
859  b = dtipsq*dtisq*w
860  IF( c.EQ.zero ) THEN
861  IF( a.EQ.zero ) THEN
862  IF( .NOT.swtch ) THEN
863  IF( orgati ) THEN
864  a = z( i )*z( i ) + dtipsq*dtipsq*
865  $ ( dpsi+dphi )
866  ELSE
867  a = z( ip1 )*z( ip1 ) +
868  $ dtisq*dtisq*( dpsi+dphi )
869  END IF
870  ELSE
871  a = dtisq*dtisq*dpsi + dtipsq*dtipsq*dphi
872  END IF
873  END IF
874  eta = b / a
875  ELSE IF( a.LE.zero ) THEN
876  eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
877  ELSE
878  eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
879  END IF
880  ELSE
881 *
882 * Interpolation using THREE most relevant poles
883 *
884  dtiim = work( iim1 )*delta( iim1 )
885  dtiip = work( iip1 )*delta( iip1 )
886  temp = rhoinv + psi + phi
887  IF( swtch ) THEN
888  c = temp - dtiim*dpsi - dtiip*dphi
889  zz( 1 ) = dtiim*dtiim*dpsi
890  zz( 3 ) = dtiip*dtiip*dphi
891  ELSE
892  IF( orgati ) THEN
893  temp1 = z( iim1 ) / dtiim
894  temp1 = temp1*temp1
895  temp2 = ( d( iim1 )-d( iip1 ) )*
896  $ ( d( iim1 )+d( iip1 ) )*temp1
897  c = temp - dtiip*( dpsi+dphi ) - temp2
898  zz( 1 ) = z( iim1 )*z( iim1 )
899  IF( dpsi.LT.temp1 ) THEN
900  zz( 3 ) = dtiip*dtiip*dphi
901  ELSE
902  zz( 3 ) = dtiip*dtiip*( ( dpsi-temp1 )+dphi )
903  END IF
904  ELSE
905  temp1 = z( iip1 ) / dtiip
906  temp1 = temp1*temp1
907  temp2 = ( d( iip1 )-d( iim1 ) )*
908  $ ( d( iim1 )+d( iip1 ) )*temp1
909  c = temp - dtiim*( dpsi+dphi ) - temp2
910  IF( dphi.LT.temp1 ) THEN
911  zz( 1 ) = dtiim*dtiim*dpsi
912  ELSE
913  zz( 1 ) = dtiim*dtiim*( dpsi+( dphi-temp1 ) )
914  END IF
915  zz( 3 ) = z( iip1 )*z( iip1 )
916  END IF
917  END IF
918  dd( 1 ) = dtiim
919  dd( 2 ) = delta( ii )*work( ii )
920  dd( 3 ) = dtiip
921  CALL dlaed6( niter, orgati, c, dd, zz, w, eta, info )
922 *
923  IF( info.NE.0 ) THEN
924 *
925 * If INFO is not 0, i.e., DLAED6 failed, switch
926 * back to two pole interpolation
927 *
928  swtch3 = .false.
929  info = 0
930  dtipsq = work( ip1 )*delta( ip1 )
931  dtisq = work( i )*delta( i )
932  IF( .NOT.swtch ) THEN
933  IF( orgati ) THEN
934  c = w - dtipsq*dw + delsq*( z( i )/dtisq )**2
935  ELSE
936  c = w - dtisq*dw - delsq*( z( ip1 )/dtipsq )**2
937  END IF
938  ELSE
939  temp = z( ii ) / ( work( ii )*delta( ii ) )
940  IF( orgati ) THEN
941  dpsi = dpsi + temp*temp
942  ELSE
943  dphi = dphi + temp*temp
944  END IF
945  c = w - dtisq*dpsi - dtipsq*dphi
946  END IF
947  a = ( dtipsq+dtisq )*w - dtipsq*dtisq*dw
948  b = dtipsq*dtisq*w
949  IF( c.EQ.zero ) THEN
950  IF( a.EQ.zero ) THEN
951  IF( .NOT.swtch ) THEN
952  IF( orgati ) THEN
953  a = z( i )*z( i ) + dtipsq*dtipsq*
954  $ ( dpsi+dphi )
955  ELSE
956  a = z( ip1 )*z( ip1 ) +
957  $ dtisq*dtisq*( dpsi+dphi )
958  END IF
959  ELSE
960  a = dtisq*dtisq*dpsi + dtipsq*dtipsq*dphi
961  END IF
962  END IF
963  eta = b / a
964  ELSE IF( a.LE.zero ) THEN
965  eta = ( a-sqrt( abs( a*a-four*b*c ) ) ) / ( two*c )
966  ELSE
967  eta = two*b / ( a+sqrt( abs( a*a-four*b*c ) ) )
968  END IF
969  END IF
970  END IF
971 *
972 * Note, eta should be positive if w is negative, and
973 * eta should be negative otherwise. However,
974 * if for some reason caused by roundoff, eta*w > 0,
975 * we simply use one Newton step instead. This way
976 * will guarantee eta*w < 0.
977 *
978  IF( w*eta.GE.zero )
979  $ eta = -w / dw
980 *
981  eta = eta / ( sigma+sqrt( sigma*sigma+eta ) )
982  temp=tau+eta
983  IF( temp.GT.sgub .OR. temp.LT.sglb ) THEN
984  IF( w.LT.zero ) THEN
985  eta = ( sgub-tau ) / two
986  ELSE
987  eta = ( sglb-tau ) / two
988  END IF
989  IF( geomavg ) THEN
990  IF( w .LT. zero ) THEN
991  IF( tau .GT. zero ) THEN
992  eta = sqrt(sgub*tau)-tau
993  END IF
994  ELSE
995  IF( sglb .GT. zero ) THEN
996  eta = sqrt(sglb*tau)-tau
997  END IF
998  END IF
999  END IF
1000  END IF
1001 *
1002  prew = w
1003 *
1004  tau = tau + eta
1005  sigma = sigma + eta
1006 *
1007  DO 200 j = 1, n
1008  work( j ) = work( j ) + eta
1009  delta( j ) = delta( j ) - eta
1010  200 CONTINUE
1011 *
1012 * Evaluate PSI and the derivative DPSI
1013 *
1014  dpsi = zero
1015  psi = zero
1016  erretm = zero
1017  DO 210 j = 1, iim1
1018  temp = z( j ) / ( work( j )*delta( j ) )
1019  psi = psi + z( j )*temp
1020  dpsi = dpsi + temp*temp
1021  erretm = erretm + psi
1022  210 CONTINUE
1023  erretm = abs( erretm )
1024 *
1025 * Evaluate PHI and the derivative DPHI
1026 *
1027  dphi = zero
1028  phi = zero
1029  DO 220 j = n, iip1, -1
1030  temp = z( j ) / ( work( j )*delta( j ) )
1031  phi = phi + z( j )*temp
1032  dphi = dphi + temp*temp
1033  erretm = erretm + phi
1034  220 CONTINUE
1035 *
1036  tau2 = work( ii )*delta( ii )
1037  temp = z( ii ) / tau2
1038  dw = dpsi + dphi + temp*temp
1039  temp = z( ii )*temp
1040  w = rhoinv + phi + psi + temp
1041  erretm = eight*( phi-psi ) + erretm + two*rhoinv
1042  $ + three*abs( temp )
1043 * $ + ABS( TAU2 )*DW
1044 *
1045  IF( w*prew.GT.zero .AND. abs( w ).GT.abs( prew ) / ten )
1046  $ swtch = .NOT.swtch
1047 *
1048  230 CONTINUE
1049 *
1050 * Return with INFO = 1, NITER = MAXIT and not converged
1051 *
1052  info = 1
1053 *
1054  END IF
1055 *
1056  240 CONTINUE
1057  RETURN
1058 *
1059 * End of DLASD4
1060 *
1061  END
dlaed6
subroutine dlaed6(KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)
DLAED6 used by sstedc. Computes one Newton step in solution of the secular equation.
Definition: dlaed6.f:142
dlasd5
subroutine dlasd5(I, D, Z, DELTA, RHO, DSIGMA, WORK)
DLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification ...
Definition: dlasd5.f:118
dlasd4
subroutine dlasd4(N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO)
DLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modif...
Definition: dlasd4.f:155