LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlasd4()

subroutine dlasd4 ( integer n,
integer i,
double precision, dimension( * ) d,
double precision, dimension( * ) z,
double precision, dimension( * ) delta,
double precision rho,
double precision sigma,
double precision, dimension( * ) work,
integer info )

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.

Download DLASD4 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> This subroutine computes the square root of the I-th updated
!> eigenvalue of a positive symmetric rank-one modification to
!> a positive diagonal matrix whose entries are given as the squares
!> of the corresponding entries in the array d, and that
!>
!>        0 <= D(i) < D(j)  for  i < j
!>
!> and that RHO > 0. This is arranged by the calling routine, and is
!> no loss in generality.  The rank-one modified system is thus
!>
!>        diag( D ) * diag( D ) +  RHO * Z * Z_transpose.
!>
!> where we assume the Euclidean norm of Z is 1.
!>
!> The method consists of approximating the rational functions in the
!> secular equation by simpler interpolating rational functions.
!>
Parameters
[in]N
!>          N is INTEGER
!>         The length of all arrays.
!>
[in]I
!>          I is INTEGER
!>         The index of the eigenvalue to be computed.  1 <= I <= N.
!>
[in]D
!>          D is DOUBLE PRECISION array, dimension ( N )
!>         The original eigenvalues.  It is assumed that they are in
!>         order, 0 <= D(I) < D(J)  for I < J.
!>
[in]Z
!>          Z is DOUBLE PRECISION array, dimension ( N )
!>         The components of the updating vector.
!>
[out]DELTA
!>          DELTA is DOUBLE PRECISION array, dimension ( N )
!>         If N .ne. 1, DELTA contains (D(j) - sigma_I) in its  j-th
!>         component.  If N = 1, then DELTA(1) = 1.  The vector DELTA
!>         contains the information necessary to construct the
!>         (singular) eigenvectors.
!>
[in]RHO
!>          RHO is DOUBLE PRECISION
!>         The scalar in the symmetric updating formula.
!>
[out]SIGMA
!>          SIGMA is DOUBLE PRECISION
!>         The computed sigma_I, the I-th updated eigenvalue.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension ( N )
!>         If N .ne. 1, WORK contains (D(j) + sigma_I) in its  j-th
!>         component.  If N = 1, then WORK( 1 ) = 1.
!>
[out]INFO
!>          INFO is INTEGER
!>         = 0:  successful exit
!>         > 0:  if INFO = 1, the updating process failed.
!>
Internal Parameters:
!>  Logical variable ORGATI (origin-at-i?) is used for distinguishing
!>  whether D(i) or D(i+1) is treated as the origin.
!>
!>            ORGATI = .true.    origin at i
!>            ORGATI = .false.   origin at i+1
!>
!>  Logical variable SWTCH3 (switch-for-3-poles?) is for noting
!>  if we are working with THREE poles!
!>
!>  MAXIT is the maximum number of iterations allowed for each
!>  eigenvalue.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

Definition at line 150 of file dlasd4.f.

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 DOUBLE PRECISION RHO, SIGMA
159* ..
160* .. Array Arguments ..
161 DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), Z( * )
162* ..
163*
164* =====================================================================
165*
166* .. Parameters ..
167 INTEGER MAXIT
168 parameter( maxit = 400 )
169 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
170 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
171 $ three = 3.0d+0, four = 4.0d+0, eight = 8.0d+0,
172 $ ten = 10.0d+0 )
173* ..
174* .. Local Scalars ..
175 LOGICAL ORGATI, SWTCH, SWTCH3, GEOMAVG
176 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
177 DOUBLE PRECISION 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 DOUBLE PRECISION DD( 3 ), ZZ( 3 )
184* ..
185* .. External Subroutines ..
186 EXTERNAL dlaed6, dlasd5
187* ..
188* .. External Functions ..
189 DOUBLE PRECISION DLAMCH
190 EXTERNAL dlamch
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 dlasd5( i, d, z, delta, rho, sigma, work )
214 RETURN
215 END IF
216*
217* Compute machine epsilon
218*
219 eps = dlamch( '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 dlaed6( niter, orgati, c, dd, zz, w, eta, info )
698*
699 IF( info.NE.0 ) THEN
700*
701* If INFO is not 0, i.e., DLAED6 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 dlaed6( niter, orgati, c, dd, zz, w, eta, info )
917*
918 IF( info.NE.0 ) THEN
919*
920* If INFO is not 0, i.e., DLAED6 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 DLASD4
1055*
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:139
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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:114
Here is the call graph for this function:
Here is the caller graph for this function: