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