LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slasd4()

subroutine slasd4 ( integer  N,
integer  I,
real, dimension( * )  D,
real, dimension( * )  Z,
real, dimension( * )  DELTA,
real  RHO,
real  SIGMA,
real, dimension( * )  WORK,
integer  INFO 
)

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.

Download SLASD4 + 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 REAL 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 REAL array, dimension ( N )
         The components of the updating vector.
[out]DELTA
          DELTA is REAL 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 REAL
         The scalar in the symmetric updating formula.
[out]SIGMA
          SIGMA is REAL
         The computed sigma_I, the I-th updated eigenvalue.
[out]WORK
          WORK is REAL 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 152 of file slasd4.f.

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 *
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
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: