LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \ingroup lasd4
144*
145*> \par Contributors:
146* ==================
147*>
148*> Ren-Cang Li, Computer Science Division, University of California
149*> at Berkeley, USA
150*>
151* =====================================================================
152 SUBROUTINE dlasd4( 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 DOUBLE PRECISION RHO, SIGMA
161* ..
162* .. Array Arguments ..
163 DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), Z( * )
164* ..
165*
166* =====================================================================
167*
168* .. Parameters ..
169 INTEGER MAXIT
170 parameter( maxit = 400 )
171 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
172 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
173 $ three = 3.0d+0, four = 4.0d+0, eight = 8.0d+0,
174 $ ten = 10.0d+0 )
175* ..
176* .. Local Scalars ..
177 LOGICAL ORGATI, SWTCH, SWTCH3, GEOMAVG
178 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
179 DOUBLE PRECISION 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 DOUBLE PRECISION DD( 3 ), ZZ( 3 )
186* ..
187* .. External Subroutines ..
188 EXTERNAL dlaed6, dlasd5
189* ..
190* .. External Functions ..
191 DOUBLE PRECISION DLAMCH
192 EXTERNAL dlamch
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 dlasd5( i, d, z, delta, rho, sigma, work )
216 RETURN
217 END IF
218*
219* Compute machine epsilon
220*
221 eps = dlamch( '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 dlaed6( niter, orgati, c, dd, zz, w, eta, info )
700*
701 IF( info.NE.0 ) THEN
702*
703* If INFO is not 0, i.e., DLAED6 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 dlaed6( niter, orgati, c, dd, zz, w, eta, info )
919*
920 IF( info.NE.0 ) THEN
921*
922* If INFO is not 0, i.e., DLAED6 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 DLASD4
1057*
1058 END
dlaed6
subroutine dlaed6(kniter, orgati, rho, d, z, finit, tau, info)
DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation.
Definition dlaed6.f:140
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:153
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:116