OpenRadioss 2025.1.11
OpenRadioss project
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*> \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
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:116
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:153
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:140
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21