LAPACK  3.6.1
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 *> \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 slasd4( 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  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  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 slaed6( niter, orgati, c, dd, zz, w, eta, info )
703 *
704  IF( info.NE.0 ) THEN
705 *
706 * If INFO is not 0, i.e., SLAED6 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 slaed6( niter, orgati, c, dd, zz, w, eta, info )
922 *
923  IF( info.NE.0 ) THEN
924 *
925 * If INFO is not 0, i.e., SLAED6 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 SLASD4
1060 *
1061  END
slasd4
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:155
slaed6
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:142
slasd5
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:118