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