1 SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER I, INFO, N
10 DOUBLE PRECISION RHO, SIGMA
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), Z( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * This subroutine computes the square root of the I-th updated
20 * eigenvalue of a positive symmetric rank-one modification to
21 * a positive diagonal matrix whose entries are given as the squares
22 * of the corresponding entries in the array d, and that
23 *
24 * 0 <= D(i) < D(j) for i < j
25 *
26 * and that RHO > 0. This is arranged by the calling routine, and is
27 * no loss in generality. The rank-one modified system is thus
28 *
29 * diag( D ) * diag( D ) + RHO * Z * Z_transpose.
30 *
31 * where we assume the Euclidean norm of Z is 1.
32 *
33 * The method consists of approximating the rational functions in the
34 * secular equation by simpler interpolating rational functions.
35 *
36 * Arguments
37 * =========
38 *
39 * N (input) INTEGER
40 * The length of all arrays.
41 *
42 * I (input) INTEGER
43 * The index of the eigenvalue to be computed. 1 <= I <= N.
44 *
45 * D (input) DOUBLE PRECISION array, dimension ( N )
46 * The original eigenvalues. It is assumed that they are in
47 * order, 0 <= D(I) < D(J) for I < J.
48 *
49 * Z (input) DOUBLE PRECISION array, dimension ( N )
50 * The components of the updating vector.
51 *
52 * DELTA (output) DOUBLE PRECISION array, dimension ( N )
53 * If N .ne. 1, DELTA contains (D(j) - sigma_I) in its j-th
54 * component. If N = 1, then DELTA(1) = 1. The vector DELTA
55 * contains the information necessary to construct the
56 * (singular) eigenvectors.
57 *
58 * RHO (input) DOUBLE PRECISION
59 * The scalar in the symmetric updating formula.
60 *
61 * SIGMA (output) DOUBLE PRECISION
62 * The computed sigma_I, the I-th updated eigenvalue.
63 *
64 * WORK (workspace) DOUBLE PRECISION array, dimension ( N )
65 * If N .ne. 1, WORK contains (D(j) + sigma_I) in its j-th
66 * component. If N = 1, then WORK( 1 ) = 1.
67 *
68 * INFO (output) INTEGER
69 * = 0: successful exit
70 * > 0: if INFO = 1, the updating process failed.
71 *
72 * Internal Parameters
73 * ===================
74 *
75 * Logical variable ORGATI (origin-at-i?) is used for distinguishing
76 * whether D(i) or D(i+1) is treated as the origin.
77 *
78 * ORGATI = .true. origin at i
79 * ORGATI = .false. origin at i+1
80 *
81 * Logical variable SWTCH3 (switch-for-3-poles?) is for noting
82 * if we are working with THREE poles!
83 *
84 * MAXIT is the maximum number of iterations allowed for each
85 * eigenvalue.
86 *
87 * Further Details
88 * ===============
89 *
90 * Based on contributions by
91 * Ren-Cang Li, Computer Science Division, University of California
92 * at Berkeley, USA
93 *
94 * =====================================================================
95 *
96 * .. Parameters ..
97 INTEGER MAXIT
98 PARAMETER ( MAXIT = 64 )
99 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
100 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
101 $ THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0,
102 $ TEN = 10.0D+0 )
103 * ..
104 * .. Local Scalars ..
105 LOGICAL ORGATI, SWTCH, SWTCH3
106 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
107 DOUBLE PRECISION A, B, C, DELSQ, DELSQ2, DPHI, DPSI, DTIIM,
108 $ DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS,
109 $ ERRETM, ETA, PHI, PREW, PSI, RHOINV, SG2LB,
110 $ SG2UB, TAU, TEMP, TEMP1, TEMP2, W
111 * ..
112 * .. Local Arrays ..
113 DOUBLE PRECISION DD( 3 ), ZZ( 3 )
114 * ..
115 * .. External Subroutines ..
116 EXTERNAL DLAED6, DLASD5
117 * ..
118 * .. External Functions ..
119 DOUBLE PRECISION DLAMCH
120 EXTERNAL DLAMCH
121 * ..
122 * .. Intrinsic Functions ..
123 INTRINSIC ABS, MAX, MIN, SQRT
124 * ..
125 * .. Executable Statements ..
126 *
127 * Since this routine is called in an inner loop, we do no argument
128 * checking.
129 *
130 * Quick return for N=1 and 2.
131 *
132 INFO = 0
133 IF( N.EQ.1 ) THEN
134 *
135 * Presumably, I=1 upon entry
136 *
137 SIGMA = SQRT( D( 1 )*D( 1 )+RHO*Z( 1 )*Z( 1 ) )
138 DELTA( 1 ) = ONE
139 WORK( 1 ) = ONE
140 RETURN
141 END IF
142 IF( N.EQ.2 ) THEN
143 CALL DLASD5( I, D, Z, DELTA, RHO, SIGMA, WORK )
144 RETURN
145 END IF
146 *
147 * Compute machine epsilon
148 *
149 EPS = DLAMCH( 'Epsilon' )
150 RHOINV = ONE / RHO
151 *
152 * The case I = N
153 *
154 IF( I.EQ.N ) THEN
155 *
156 * Initialize some basic variables
157 *
158 II = N - 1
159 NITER = 1
160 *
161 * Calculate initial guess
162 *
163 TEMP = RHO / TWO
164 *
165 * If ||Z||_2 is not one, then TEMP should be set to
166 * RHO * ||Z||_2^2 / TWO
167 *
168 TEMP1 = TEMP / ( D( N )+SQRT( D( N )*D( N )+TEMP ) )
169 DO 10 J = 1, N
170 WORK( J ) = D( J ) + D( N ) + TEMP1
171 DELTA( J ) = ( D( J )-D( N ) ) - TEMP1
172 10 CONTINUE
173 *
174 PSI = ZERO
175 DO 20 J = 1, N - 2
176 PSI = PSI + Z( J )*Z( J ) / ( DELTA( J )*WORK( J ) )
177 20 CONTINUE
178 *
179 C = RHOINV + PSI
180 W = C + Z( II )*Z( II ) / ( DELTA( II )*WORK( II ) ) +
181 $ Z( N )*Z( N ) / ( DELTA( N )*WORK( N ) )
182 *
183 IF( W.LE.ZERO ) THEN
184 TEMP1 = SQRT( D( N )*D( N )+RHO )
185 TEMP = Z( N-1 )*Z( N-1 ) / ( ( D( N-1 )+TEMP1 )*
186 $ ( D( N )-D( N-1 )+RHO / ( D( N )+TEMP1 ) ) ) +
187 $ Z( N )*Z( N ) / RHO
188 *
189 * The following TAU is to approximate
190 * SIGMA_n^2 - D( N )*D( N )
191 *
192 IF( C.LE.TEMP ) THEN
193 TAU = RHO
194 ELSE
195 DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
196 A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
197 B = Z( N )*Z( N )*DELSQ
198 IF( A.LT.ZERO ) THEN
199 TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
200 ELSE
201 TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
202 END IF
203 END IF
204 *
205 * It can be proved that
206 * D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU <= D(N)^2+RHO
207 *
208 ELSE
209 DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
210 A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
211 B = Z( N )*Z( N )*DELSQ
212 *
213 * The following TAU is to approximate
214 * SIGMA_n^2 - D( N )*D( N )
215 *
216 IF( A.LT.ZERO ) THEN
217 TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
218 ELSE
219 TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
220 END IF
221 *
222 * It can be proved that
223 * D(N)^2 < D(N)^2+TAU < SIGMA(N)^2 < D(N)^2+RHO/2
224 *
225 END IF
226 *
227 * The following ETA is to approximate SIGMA_n - D( N )
228 *
229 ETA = TAU / ( D( N )+SQRT( D( N )*D( N )+TAU ) )
230 *
231 SIGMA = D( N ) + ETA
232 DO 30 J = 1, N
233 DELTA( J ) = ( D( J )-D( I ) ) - ETA
234 WORK( J ) = D( J ) + D( I ) + ETA
235 30 CONTINUE
236 *
237 * Evaluate PSI and the derivative DPSI
238 *
239 DPSI = ZERO
240 PSI = ZERO
241 ERRETM = ZERO
242 DO 40 J = 1, II
243 TEMP = Z( J ) / ( DELTA( J )*WORK( J ) )
244 PSI = PSI + Z( J )*TEMP
245 DPSI = DPSI + TEMP*TEMP
246 ERRETM = ERRETM + PSI
247 40 CONTINUE
248 ERRETM = ABS( ERRETM )
249 *
250 * Evaluate PHI and the derivative DPHI
251 *
252 TEMP = Z( N ) / ( DELTA( N )*WORK( N ) )
253 PHI = Z( N )*TEMP
254 DPHI = TEMP*TEMP
255 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
256 $ ABS( TAU )*( DPSI+DPHI )
257 *
258 W = RHOINV + PHI + PSI
259 *
260 * Test for convergence
261 *
262 IF( ABS( W ).LE.EPS*ERRETM ) THEN
263 GO TO 240
264 END IF
265 *
266 * Calculate the new step
267 *
268 NITER = NITER + 1
269 DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
270 DTNSQ = WORK( N )*DELTA( N )
271 C = W - DTNSQ1*DPSI - DTNSQ*DPHI
272 A = ( DTNSQ+DTNSQ1 )*W - DTNSQ*DTNSQ1*( DPSI+DPHI )
273 B = DTNSQ*DTNSQ1*W
274 IF( C.LT.ZERO )
275 $ C = ABS( C )
276 IF( C.EQ.ZERO ) THEN
277 ETA = RHO - SIGMA*SIGMA
278 ELSE IF( A.GE.ZERO ) THEN
279 ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
280 ELSE
281 ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
282 END IF
283 *
284 * Note, eta should be positive if w is negative, and
285 * eta should be negative otherwise. However,
286 * if for some reason caused by roundoff, eta*w > 0,
287 * we simply use one Newton step instead. This way
288 * will guarantee eta*w < 0.
289 *
290 IF( W*ETA.GT.ZERO )
291 $ ETA = -W / ( DPSI+DPHI )
292 TEMP = ETA - DTNSQ
293 IF( TEMP.GT.RHO )
294 $ ETA = RHO + DTNSQ
295 *
296 TAU = TAU + ETA
297 ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
298 DO 50 J = 1, N
299 DELTA( J ) = DELTA( J ) - ETA
300 WORK( J ) = WORK( J ) + ETA
301 50 CONTINUE
302 *
303 SIGMA = SIGMA + ETA
304 *
305 * Evaluate PSI and the derivative DPSI
306 *
307 DPSI = ZERO
308 PSI = ZERO
309 ERRETM = ZERO
310 DO 60 J = 1, II
311 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
312 PSI = PSI + Z( J )*TEMP
313 DPSI = DPSI + TEMP*TEMP
314 ERRETM = ERRETM + PSI
315 60 CONTINUE
316 ERRETM = ABS( ERRETM )
317 *
318 * Evaluate PHI and the derivative DPHI
319 *
320 TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )
321 PHI = Z( N )*TEMP
322 DPHI = TEMP*TEMP
323 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
324 $ ABS( TAU )*( DPSI+DPHI )
325 *
326 W = RHOINV + PHI + PSI
327 *
328 * Main loop to update the values of the array DELTA
329 *
330 ITER = NITER + 1
331 *
332 DO 90 NITER = ITER, MAXIT
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 DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
343 DTNSQ = WORK( N )*DELTA( N )
344 C = W - DTNSQ1*DPSI - DTNSQ*DPHI
345 A = ( DTNSQ+DTNSQ1 )*W - DTNSQ1*DTNSQ*( DPSI+DPHI )
346 B = DTNSQ1*DTNSQ*W
347 IF( A.GE.ZERO ) THEN
348 ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
349 ELSE
350 ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
351 END IF
352 *
353 * Note, eta should be positive if w is negative, and
354 * eta should be negative otherwise. However,
355 * if for some reason caused by roundoff, eta*w > 0,
356 * we simply use one Newton step instead. This way
357 * will guarantee eta*w < 0.
358 *
359 IF( W*ETA.GT.ZERO )
360 $ ETA = -W / ( DPSI+DPHI )
361 TEMP = ETA - DTNSQ
362 IF( TEMP.LE.ZERO )
363 $ ETA = ETA / TWO
364 *
365 TAU = TAU + ETA
366 ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
367 DO 70 J = 1, N
368 DELTA( J ) = DELTA( J ) - ETA
369 WORK( J ) = WORK( J ) + ETA
370 70 CONTINUE
371 *
372 SIGMA = SIGMA + ETA
373 *
374 * Evaluate PSI and the derivative DPSI
375 *
376 DPSI = ZERO
377 PSI = ZERO
378 ERRETM = ZERO
379 DO 80 J = 1, II
380 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
381 PSI = PSI + Z( J )*TEMP
382 DPSI = DPSI + TEMP*TEMP
383 ERRETM = ERRETM + PSI
384 80 CONTINUE
385 ERRETM = ABS( ERRETM )
386 *
387 * Evaluate PHI and the derivative DPHI
388 *
389 TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )
390 PHI = Z( N )*TEMP
391 DPHI = TEMP*TEMP
392 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
393 $ ABS( TAU )*( DPSI+DPHI )
394 *
395 W = RHOINV + PHI + PSI
396 90 CONTINUE
397 *
398 * Return with INFO = 1, NITER = MAXIT and not converged
399 *
400 INFO = 1
401 GO TO 240
402 *
403 * End for the case I = N
404 *
405 ELSE
406 *
407 * The case for I < N
408 *
409 NITER = 1
410 IP1 = I + 1
411 *
412 * Calculate initial guess
413 *
414 DELSQ = ( D( IP1 )-D( I ) )*( D( IP1 )+D( I ) )
415 DELSQ2 = DELSQ / TWO
416 TEMP = DELSQ2 / ( D( I )+SQRT( D( I )*D( I )+DELSQ2 ) )
417 DO 100 J = 1, N
418 WORK( J ) = D( J ) + D( I ) + TEMP
419 DELTA( J ) = ( D( J )-D( I ) ) - TEMP
420 100 CONTINUE
421 *
422 PSI = ZERO
423 DO 110 J = 1, I - 1
424 PSI = PSI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
425 110 CONTINUE
426 *
427 PHI = ZERO
428 DO 120 J = N, I + 2, -1
429 PHI = PHI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
430 120 CONTINUE
431 C = RHOINV + PSI + PHI
432 W = C + Z( I )*Z( I ) / ( WORK( I )*DELTA( I ) ) +
433 $ Z( IP1 )*Z( IP1 ) / ( WORK( IP1 )*DELTA( IP1 ) )
434 *
435 IF( W.GT.ZERO ) THEN
436 *
437 * d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2
438 *
439 * We choose d(i) as origin.
440 *
441 ORGATI = .TRUE.
442 SG2LB = ZERO
443 SG2UB = DELSQ2
444 A = C*DELSQ + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
445 B = Z( I )*Z( I )*DELSQ
446 IF( A.GT.ZERO ) THEN
447 TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
448 ELSE
449 TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
450 END IF
451 *
452 * TAU now is an estimation of SIGMA^2 - D( I )^2. The
453 * following, however, is the corresponding estimation of
454 * SIGMA - D( I ).
455 *
456 ETA = TAU / ( D( I )+SQRT( D( I )*D( I )+TAU ) )
457 ELSE
458 *
459 * (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2
460 *
461 * We choose d(i+1) as origin.
462 *
463 ORGATI = .FALSE.
464 SG2LB = -DELSQ2
465 SG2UB = ZERO
466 A = C*DELSQ - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
467 B = Z( IP1 )*Z( IP1 )*DELSQ
468 IF( A.LT.ZERO ) THEN
469 TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
470 ELSE
471 TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
472 END IF
473 *
474 * TAU now is an estimation of SIGMA^2 - D( IP1 )^2. The
475 * following, however, is the corresponding estimation of
476 * SIGMA - D( IP1 ).
477 *
478 ETA = TAU / ( D( IP1 )+SQRT( ABS( D( IP1 )*D( IP1 )+
479 $ TAU ) ) )
480 END IF
481 *
482 IF( ORGATI ) THEN
483 II = I
484 SIGMA = D( I ) + ETA
485 DO 130 J = 1, N
486 WORK( J ) = D( J ) + D( I ) + ETA
487 DELTA( J ) = ( D( J )-D( I ) ) - ETA
488 130 CONTINUE
489 ELSE
490 II = I + 1
491 SIGMA = D( IP1 ) + ETA
492 DO 140 J = 1, N
493 WORK( J ) = D( J ) + D( IP1 ) + ETA
494 DELTA( J ) = ( D( J )-D( IP1 ) ) - ETA
495 140 CONTINUE
496 END IF
497 IIM1 = II - 1
498 IIP1 = II + 1
499 *
500 * Evaluate PSI and the derivative DPSI
501 *
502 DPSI = ZERO
503 PSI = ZERO
504 ERRETM = ZERO
505 DO 150 J = 1, IIM1
506 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
507 PSI = PSI + Z( J )*TEMP
508 DPSI = DPSI + TEMP*TEMP
509 ERRETM = ERRETM + PSI
510 150 CONTINUE
511 ERRETM = ABS( ERRETM )
512 *
513 * Evaluate PHI and the derivative DPHI
514 *
515 DPHI = ZERO
516 PHI = ZERO
517 DO 160 J = N, IIP1, -1
518 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
519 PHI = PHI + Z( J )*TEMP
520 DPHI = DPHI + TEMP*TEMP
521 ERRETM = ERRETM + PHI
522 160 CONTINUE
523 *
524 W = RHOINV + PHI + PSI
525 *
526 * W is the value of the secular function with
527 * its ii-th element removed.
528 *
529 SWTCH3 = .FALSE.
530 IF( ORGATI ) THEN
531 IF( W.LT.ZERO )
532 $ SWTCH3 = .TRUE.
533 ELSE
534 IF( W.GT.ZERO )
535 $ SWTCH3 = .TRUE.
536 END IF
537 IF( II.EQ.1 .OR. II.EQ.N )
538 $ SWTCH3 = .FALSE.
539 *
540 TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
541 DW = DPSI + DPHI + TEMP*TEMP
542 TEMP = Z( II )*TEMP
543 W = W + TEMP
544 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
545 $ THREE*ABS( TEMP ) + ABS( TAU )*DW
546 *
547 * Test for convergence
548 *
549 IF( ABS( W ).LE.EPS*ERRETM ) THEN
550 GO TO 240
551 END IF
552 *
553 IF( W.LE.ZERO ) THEN
554 SG2LB = MAX( SG2LB, TAU )
555 ELSE
556 SG2UB = MIN( SG2UB, TAU )
557 END IF
558 *
559 * Calculate the new step
560 *
561 NITER = NITER + 1
562 IF( .NOT.SWTCH3 ) THEN
563 DTIPSQ = WORK( IP1 )*DELTA( IP1 )
564 DTISQ = WORK( I )*DELTA( I )
565 IF( ORGATI ) THEN
566 C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
567 ELSE
568 C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
569 END IF
570 A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
571 B = DTIPSQ*DTISQ*W
572 IF( C.EQ.ZERO ) THEN
573 IF( A.EQ.ZERO ) THEN
574 IF( ORGATI ) THEN
575 A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*( DPSI+DPHI )
576 ELSE
577 A = Z( IP1 )*Z( IP1 ) + DTISQ*DTISQ*( DPSI+DPHI )
578 END IF
579 END IF
580 ETA = B / A
581 ELSE IF( A.LE.ZERO ) THEN
582 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
583 ELSE
584 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
585 END IF
586 ELSE
587 *
588 * Interpolation using THREE most relevant poles
589 *
590 DTIIM = WORK( IIM1 )*DELTA( IIM1 )
591 DTIIP = WORK( IIP1 )*DELTA( IIP1 )
592 TEMP = RHOINV + PSI + PHI
593 IF( ORGATI ) THEN
594 TEMP1 = Z( IIM1 ) / DTIIM
595 TEMP1 = TEMP1*TEMP1
596 C = ( TEMP - DTIIP*( DPSI+DPHI ) ) -
597 $ ( D( IIM1 )-D( IIP1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
598 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
599 IF( DPSI.LT.TEMP1 ) THEN
600 ZZ( 3 ) = DTIIP*DTIIP*DPHI
601 ELSE
602 ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
603 END IF
604 ELSE
605 TEMP1 = Z( IIP1 ) / DTIIP
606 TEMP1 = TEMP1*TEMP1
607 C = ( TEMP - DTIIM*( DPSI+DPHI ) ) -
608 $ ( D( IIP1 )-D( IIM1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
609 IF( DPHI.LT.TEMP1 ) THEN
610 ZZ( 1 ) = DTIIM*DTIIM*DPSI
611 ELSE
612 ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
613 END IF
614 ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
615 END IF
616 ZZ( 2 ) = Z( II )*Z( II )
617 DD( 1 ) = DTIIM
618 DD( 2 ) = DELTA( II )*WORK( II )
619 DD( 3 ) = DTIIP
620 CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
621 IF( INFO.NE.0 )
622 $ GO TO 240
623 END IF
624 *
625 * Note, eta should be positive if w is negative, and
626 * eta should be negative otherwise. However,
627 * if for some reason caused by roundoff, eta*w > 0,
628 * we simply use one Newton step instead. This way
629 * will guarantee eta*w < 0.
630 *
631 IF( W*ETA.GE.ZERO )
632 $ ETA = -W / DW
633 IF( ORGATI ) THEN
634 TEMP1 = WORK( I )*DELTA( I )
635 TEMP = ETA - TEMP1
636 ELSE
637 TEMP1 = WORK( IP1 )*DELTA( IP1 )
638 TEMP = ETA - TEMP1
639 END IF
640 IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN
641 IF( W.LT.ZERO ) THEN
642 ETA = ( SG2UB-TAU ) / TWO
643 ELSE
644 ETA = ( SG2LB-TAU ) / TWO
645 END IF
646 END IF
647 *
648 TAU = TAU + ETA
649 ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
650 *
651 PREW = W
652 *
653 SIGMA = SIGMA + ETA
654 DO 170 J = 1, N
655 WORK( J ) = WORK( J ) + ETA
656 DELTA( J ) = DELTA( J ) - ETA
657 170 CONTINUE
658 *
659 * Evaluate PSI and the derivative DPSI
660 *
661 DPSI = ZERO
662 PSI = ZERO
663 ERRETM = ZERO
664 DO 180 J = 1, IIM1
665 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
666 PSI = PSI + Z( J )*TEMP
667 DPSI = DPSI + TEMP*TEMP
668 ERRETM = ERRETM + PSI
669 180 CONTINUE
670 ERRETM = ABS( ERRETM )
671 *
672 * Evaluate PHI and the derivative DPHI
673 *
674 DPHI = ZERO
675 PHI = ZERO
676 DO 190 J = N, IIP1, -1
677 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
678 PHI = PHI + Z( J )*TEMP
679 DPHI = DPHI + TEMP*TEMP
680 ERRETM = ERRETM + PHI
681 190 CONTINUE
682 *
683 TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
684 DW = DPSI + DPHI + TEMP*TEMP
685 TEMP = Z( II )*TEMP
686 W = RHOINV + PHI + PSI + TEMP
687 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
688 $ THREE*ABS( TEMP ) + ABS( TAU )*DW
689 *
690 IF( W.LE.ZERO ) THEN
691 SG2LB = MAX( SG2LB, TAU )
692 ELSE
693 SG2UB = MIN( SG2UB, TAU )
694 END IF
695 *
696 SWTCH = .FALSE.
697 IF( ORGATI ) THEN
698 IF( -W.GT.ABS( PREW ) / TEN )
699 $ SWTCH = .TRUE.
700 ELSE
701 IF( W.GT.ABS( PREW ) / TEN )
702 $ SWTCH = .TRUE.
703 END IF
704 *
705 * Main loop to update the values of the array DELTA and WORK
706 *
707 ITER = NITER + 1
708 *
709 DO 230 NITER = ITER, MAXIT
710 *
711 * Test for convergence
712 *
713 IF( ABS( W ).LE.EPS*ERRETM ) THEN
714 GO TO 240
715 END IF
716 *
717 * Calculate the new step
718 *
719 IF( .NOT.SWTCH3 ) THEN
720 DTIPSQ = WORK( IP1 )*DELTA( IP1 )
721 DTISQ = WORK( I )*DELTA( I )
722 IF( .NOT.SWTCH ) THEN
723 IF( ORGATI ) THEN
724 C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
725 ELSE
726 C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
727 END IF
728 ELSE
729 TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
730 IF( ORGATI ) THEN
731 DPSI = DPSI + TEMP*TEMP
732 ELSE
733 DPHI = DPHI + TEMP*TEMP
734 END IF
735 C = W - DTISQ*DPSI - DTIPSQ*DPHI
736 END IF
737 A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
738 B = DTIPSQ*DTISQ*W
739 IF( C.EQ.ZERO ) THEN
740 IF( A.EQ.ZERO ) THEN
741 IF( .NOT.SWTCH ) THEN
742 IF( ORGATI ) THEN
743 A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*
744 $ ( DPSI+DPHI )
745 ELSE
746 A = Z( IP1 )*Z( IP1 ) +
747 $ DTISQ*DTISQ*( DPSI+DPHI )
748 END IF
749 ELSE
750 A = DTISQ*DTISQ*DPSI + DTIPSQ*DTIPSQ*DPHI
751 END IF
752 END IF
753 ETA = B / A
754 ELSE IF( A.LE.ZERO ) THEN
755 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
756 ELSE
757 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
758 END IF
759 ELSE
760 *
761 * Interpolation using THREE most relevant poles
762 *
763 DTIIM = WORK( IIM1 )*DELTA( IIM1 )
764 DTIIP = WORK( IIP1 )*DELTA( IIP1 )
765 TEMP = RHOINV + PSI + PHI
766 IF( SWTCH ) THEN
767 C = TEMP - DTIIM*DPSI - DTIIP*DPHI
768 ZZ( 1 ) = DTIIM*DTIIM*DPSI
769 ZZ( 3 ) = DTIIP*DTIIP*DPHI
770 ELSE
771 IF( ORGATI ) THEN
772 TEMP1 = Z( IIM1 ) / DTIIM
773 TEMP1 = TEMP1*TEMP1
774 TEMP2 = ( D( IIM1 )-D( IIP1 ) )*
775 $ ( D( IIM1 )+D( IIP1 ) )*TEMP1
776 C = TEMP - DTIIP*( DPSI+DPHI ) - TEMP2
777 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
778 IF( DPSI.LT.TEMP1 ) THEN
779 ZZ( 3 ) = DTIIP*DTIIP*DPHI
780 ELSE
781 ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
782 END IF
783 ELSE
784 TEMP1 = Z( IIP1 ) / DTIIP
785 TEMP1 = TEMP1*TEMP1
786 TEMP2 = ( D( IIP1 )-D( IIM1 ) )*
787 $ ( D( IIM1 )+D( IIP1 ) )*TEMP1
788 C = TEMP - DTIIM*( DPSI+DPHI ) - TEMP2
789 IF( DPHI.LT.TEMP1 ) THEN
790 ZZ( 1 ) = DTIIM*DTIIM*DPSI
791 ELSE
792 ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
793 END IF
794 ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
795 END IF
796 END IF
797 DD( 1 ) = DTIIM
798 DD( 2 ) = DELTA( II )*WORK( II )
799 DD( 3 ) = DTIIP
800 CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
801 IF( INFO.NE.0 )
802 $ GO TO 240
803 END IF
804 *
805 * Note, eta should be positive if w is negative, and
806 * eta should be negative otherwise. However,
807 * if for some reason caused by roundoff, eta*w > 0,
808 * we simply use one Newton step instead. This way
809 * will guarantee eta*w < 0.
810 *
811 IF( W*ETA.GE.ZERO )
812 $ ETA = -W / DW
813 IF( ORGATI ) THEN
814 TEMP1 = WORK( I )*DELTA( I )
815 TEMP = ETA - TEMP1
816 ELSE
817 TEMP1 = WORK( IP1 )*DELTA( IP1 )
818 TEMP = ETA - TEMP1
819 END IF
820 IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN
821 IF( W.LT.ZERO ) THEN
822 ETA = ( SG2UB-TAU ) / TWO
823 ELSE
824 ETA = ( SG2LB-TAU ) / TWO
825 END IF
826 END IF
827 *
828 TAU = TAU + ETA
829 ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
830 *
831 SIGMA = SIGMA + ETA
832 DO 200 J = 1, N
833 WORK( J ) = WORK( J ) + ETA
834 DELTA( J ) = DELTA( J ) - ETA
835 200 CONTINUE
836 *
837 PREW = W
838 *
839 * Evaluate PSI and the derivative DPSI
840 *
841 DPSI = ZERO
842 PSI = ZERO
843 ERRETM = ZERO
844 DO 210 J = 1, IIM1
845 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
846 PSI = PSI + Z( J )*TEMP
847 DPSI = DPSI + TEMP*TEMP
848 ERRETM = ERRETM + PSI
849 210 CONTINUE
850 ERRETM = ABS( ERRETM )
851 *
852 * Evaluate PHI and the derivative DPHI
853 *
854 DPHI = ZERO
855 PHI = ZERO
856 DO 220 J = N, IIP1, -1
857 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
858 PHI = PHI + Z( J )*TEMP
859 DPHI = DPHI + TEMP*TEMP
860 ERRETM = ERRETM + PHI
861 220 CONTINUE
862 *
863 TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
864 DW = DPSI + DPHI + TEMP*TEMP
865 TEMP = Z( II )*TEMP
866 W = RHOINV + PHI + PSI + TEMP
867 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
868 $ THREE*ABS( TEMP ) + ABS( TAU )*DW
869 IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )
870 $ SWTCH = .NOT.SWTCH
871 *
872 IF( W.LE.ZERO ) THEN
873 SG2LB = MAX( SG2LB, TAU )
874 ELSE
875 SG2UB = MIN( SG2UB, TAU )
876 END IF
877 *
878 230 CONTINUE
879 *
880 * Return with INFO = 1, NITER = MAXIT and not converged
881 *
882 INFO = 1
883 *
884 END IF
885 *
886 240 CONTINUE
887 RETURN
888 *
889 * End of DLASD4
890 *
891 END
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER I, INFO, N
10 DOUBLE PRECISION RHO, SIGMA
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), Z( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * This subroutine computes the square root of the I-th updated
20 * eigenvalue of a positive symmetric rank-one modification to
21 * a positive diagonal matrix whose entries are given as the squares
22 * of the corresponding entries in the array d, and that
23 *
24 * 0 <= D(i) < D(j) for i < j
25 *
26 * and that RHO > 0. This is arranged by the calling routine, and is
27 * no loss in generality. The rank-one modified system is thus
28 *
29 * diag( D ) * diag( D ) + RHO * Z * Z_transpose.
30 *
31 * where we assume the Euclidean norm of Z is 1.
32 *
33 * The method consists of approximating the rational functions in the
34 * secular equation by simpler interpolating rational functions.
35 *
36 * Arguments
37 * =========
38 *
39 * N (input) INTEGER
40 * The length of all arrays.
41 *
42 * I (input) INTEGER
43 * The index of the eigenvalue to be computed. 1 <= I <= N.
44 *
45 * D (input) DOUBLE PRECISION array, dimension ( N )
46 * The original eigenvalues. It is assumed that they are in
47 * order, 0 <= D(I) < D(J) for I < J.
48 *
49 * Z (input) DOUBLE PRECISION array, dimension ( N )
50 * The components of the updating vector.
51 *
52 * DELTA (output) DOUBLE PRECISION array, dimension ( N )
53 * If N .ne. 1, DELTA contains (D(j) - sigma_I) in its j-th
54 * component. If N = 1, then DELTA(1) = 1. The vector DELTA
55 * contains the information necessary to construct the
56 * (singular) eigenvectors.
57 *
58 * RHO (input) DOUBLE PRECISION
59 * The scalar in the symmetric updating formula.
60 *
61 * SIGMA (output) DOUBLE PRECISION
62 * The computed sigma_I, the I-th updated eigenvalue.
63 *
64 * WORK (workspace) DOUBLE PRECISION array, dimension ( N )
65 * If N .ne. 1, WORK contains (D(j) + sigma_I) in its j-th
66 * component. If N = 1, then WORK( 1 ) = 1.
67 *
68 * INFO (output) INTEGER
69 * = 0: successful exit
70 * > 0: if INFO = 1, the updating process failed.
71 *
72 * Internal Parameters
73 * ===================
74 *
75 * Logical variable ORGATI (origin-at-i?) is used for distinguishing
76 * whether D(i) or D(i+1) is treated as the origin.
77 *
78 * ORGATI = .true. origin at i
79 * ORGATI = .false. origin at i+1
80 *
81 * Logical variable SWTCH3 (switch-for-3-poles?) is for noting
82 * if we are working with THREE poles!
83 *
84 * MAXIT is the maximum number of iterations allowed for each
85 * eigenvalue.
86 *
87 * Further Details
88 * ===============
89 *
90 * Based on contributions by
91 * Ren-Cang Li, Computer Science Division, University of California
92 * at Berkeley, USA
93 *
94 * =====================================================================
95 *
96 * .. Parameters ..
97 INTEGER MAXIT
98 PARAMETER ( MAXIT = 64 )
99 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
100 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
101 $ THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0,
102 $ TEN = 10.0D+0 )
103 * ..
104 * .. Local Scalars ..
105 LOGICAL ORGATI, SWTCH, SWTCH3
106 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
107 DOUBLE PRECISION A, B, C, DELSQ, DELSQ2, DPHI, DPSI, DTIIM,
108 $ DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS,
109 $ ERRETM, ETA, PHI, PREW, PSI, RHOINV, SG2LB,
110 $ SG2UB, TAU, TEMP, TEMP1, TEMP2, W
111 * ..
112 * .. Local Arrays ..
113 DOUBLE PRECISION DD( 3 ), ZZ( 3 )
114 * ..
115 * .. External Subroutines ..
116 EXTERNAL DLAED6, DLASD5
117 * ..
118 * .. External Functions ..
119 DOUBLE PRECISION DLAMCH
120 EXTERNAL DLAMCH
121 * ..
122 * .. Intrinsic Functions ..
123 INTRINSIC ABS, MAX, MIN, SQRT
124 * ..
125 * .. Executable Statements ..
126 *
127 * Since this routine is called in an inner loop, we do no argument
128 * checking.
129 *
130 * Quick return for N=1 and 2.
131 *
132 INFO = 0
133 IF( N.EQ.1 ) THEN
134 *
135 * Presumably, I=1 upon entry
136 *
137 SIGMA = SQRT( D( 1 )*D( 1 )+RHO*Z( 1 )*Z( 1 ) )
138 DELTA( 1 ) = ONE
139 WORK( 1 ) = ONE
140 RETURN
141 END IF
142 IF( N.EQ.2 ) THEN
143 CALL DLASD5( I, D, Z, DELTA, RHO, SIGMA, WORK )
144 RETURN
145 END IF
146 *
147 * Compute machine epsilon
148 *
149 EPS = DLAMCH( 'Epsilon' )
150 RHOINV = ONE / RHO
151 *
152 * The case I = N
153 *
154 IF( I.EQ.N ) THEN
155 *
156 * Initialize some basic variables
157 *
158 II = N - 1
159 NITER = 1
160 *
161 * Calculate initial guess
162 *
163 TEMP = RHO / TWO
164 *
165 * If ||Z||_2 is not one, then TEMP should be set to
166 * RHO * ||Z||_2^2 / TWO
167 *
168 TEMP1 = TEMP / ( D( N )+SQRT( D( N )*D( N )+TEMP ) )
169 DO 10 J = 1, N
170 WORK( J ) = D( J ) + D( N ) + TEMP1
171 DELTA( J ) = ( D( J )-D( N ) ) - TEMP1
172 10 CONTINUE
173 *
174 PSI = ZERO
175 DO 20 J = 1, N - 2
176 PSI = PSI + Z( J )*Z( J ) / ( DELTA( J )*WORK( J ) )
177 20 CONTINUE
178 *
179 C = RHOINV + PSI
180 W = C + Z( II )*Z( II ) / ( DELTA( II )*WORK( II ) ) +
181 $ Z( N )*Z( N ) / ( DELTA( N )*WORK( N ) )
182 *
183 IF( W.LE.ZERO ) THEN
184 TEMP1 = SQRT( D( N )*D( N )+RHO )
185 TEMP = Z( N-1 )*Z( N-1 ) / ( ( D( N-1 )+TEMP1 )*
186 $ ( D( N )-D( N-1 )+RHO / ( D( N )+TEMP1 ) ) ) +
187 $ Z( N )*Z( N ) / RHO
188 *
189 * The following TAU is to approximate
190 * SIGMA_n^2 - D( N )*D( N )
191 *
192 IF( C.LE.TEMP ) THEN
193 TAU = RHO
194 ELSE
195 DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
196 A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
197 B = Z( N )*Z( N )*DELSQ
198 IF( A.LT.ZERO ) THEN
199 TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
200 ELSE
201 TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
202 END IF
203 END IF
204 *
205 * It can be proved that
206 * D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU <= D(N)^2+RHO
207 *
208 ELSE
209 DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
210 A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
211 B = Z( N )*Z( N )*DELSQ
212 *
213 * The following TAU is to approximate
214 * SIGMA_n^2 - D( N )*D( N )
215 *
216 IF( A.LT.ZERO ) THEN
217 TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
218 ELSE
219 TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
220 END IF
221 *
222 * It can be proved that
223 * D(N)^2 < D(N)^2+TAU < SIGMA(N)^2 < D(N)^2+RHO/2
224 *
225 END IF
226 *
227 * The following ETA is to approximate SIGMA_n - D( N )
228 *
229 ETA = TAU / ( D( N )+SQRT( D( N )*D( N )+TAU ) )
230 *
231 SIGMA = D( N ) + ETA
232 DO 30 J = 1, N
233 DELTA( J ) = ( D( J )-D( I ) ) - ETA
234 WORK( J ) = D( J ) + D( I ) + ETA
235 30 CONTINUE
236 *
237 * Evaluate PSI and the derivative DPSI
238 *
239 DPSI = ZERO
240 PSI = ZERO
241 ERRETM = ZERO
242 DO 40 J = 1, II
243 TEMP = Z( J ) / ( DELTA( J )*WORK( J ) )
244 PSI = PSI + Z( J )*TEMP
245 DPSI = DPSI + TEMP*TEMP
246 ERRETM = ERRETM + PSI
247 40 CONTINUE
248 ERRETM = ABS( ERRETM )
249 *
250 * Evaluate PHI and the derivative DPHI
251 *
252 TEMP = Z( N ) / ( DELTA( N )*WORK( N ) )
253 PHI = Z( N )*TEMP
254 DPHI = TEMP*TEMP
255 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
256 $ ABS( TAU )*( DPSI+DPHI )
257 *
258 W = RHOINV + PHI + PSI
259 *
260 * Test for convergence
261 *
262 IF( ABS( W ).LE.EPS*ERRETM ) THEN
263 GO TO 240
264 END IF
265 *
266 * Calculate the new step
267 *
268 NITER = NITER + 1
269 DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
270 DTNSQ = WORK( N )*DELTA( N )
271 C = W - DTNSQ1*DPSI - DTNSQ*DPHI
272 A = ( DTNSQ+DTNSQ1 )*W - DTNSQ*DTNSQ1*( DPSI+DPHI )
273 B = DTNSQ*DTNSQ1*W
274 IF( C.LT.ZERO )
275 $ C = ABS( C )
276 IF( C.EQ.ZERO ) THEN
277 ETA = RHO - SIGMA*SIGMA
278 ELSE IF( A.GE.ZERO ) THEN
279 ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
280 ELSE
281 ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
282 END IF
283 *
284 * Note, eta should be positive if w is negative, and
285 * eta should be negative otherwise. However,
286 * if for some reason caused by roundoff, eta*w > 0,
287 * we simply use one Newton step instead. This way
288 * will guarantee eta*w < 0.
289 *
290 IF( W*ETA.GT.ZERO )
291 $ ETA = -W / ( DPSI+DPHI )
292 TEMP = ETA - DTNSQ
293 IF( TEMP.GT.RHO )
294 $ ETA = RHO + DTNSQ
295 *
296 TAU = TAU + ETA
297 ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
298 DO 50 J = 1, N
299 DELTA( J ) = DELTA( J ) - ETA
300 WORK( J ) = WORK( J ) + ETA
301 50 CONTINUE
302 *
303 SIGMA = SIGMA + ETA
304 *
305 * Evaluate PSI and the derivative DPSI
306 *
307 DPSI = ZERO
308 PSI = ZERO
309 ERRETM = ZERO
310 DO 60 J = 1, II
311 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
312 PSI = PSI + Z( J )*TEMP
313 DPSI = DPSI + TEMP*TEMP
314 ERRETM = ERRETM + PSI
315 60 CONTINUE
316 ERRETM = ABS( ERRETM )
317 *
318 * Evaluate PHI and the derivative DPHI
319 *
320 TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )
321 PHI = Z( N )*TEMP
322 DPHI = TEMP*TEMP
323 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
324 $ ABS( TAU )*( DPSI+DPHI )
325 *
326 W = RHOINV + PHI + PSI
327 *
328 * Main loop to update the values of the array DELTA
329 *
330 ITER = NITER + 1
331 *
332 DO 90 NITER = ITER, MAXIT
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 DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
343 DTNSQ = WORK( N )*DELTA( N )
344 C = W - DTNSQ1*DPSI - DTNSQ*DPHI
345 A = ( DTNSQ+DTNSQ1 )*W - DTNSQ1*DTNSQ*( DPSI+DPHI )
346 B = DTNSQ1*DTNSQ*W
347 IF( A.GE.ZERO ) THEN
348 ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
349 ELSE
350 ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
351 END IF
352 *
353 * Note, eta should be positive if w is negative, and
354 * eta should be negative otherwise. However,
355 * if for some reason caused by roundoff, eta*w > 0,
356 * we simply use one Newton step instead. This way
357 * will guarantee eta*w < 0.
358 *
359 IF( W*ETA.GT.ZERO )
360 $ ETA = -W / ( DPSI+DPHI )
361 TEMP = ETA - DTNSQ
362 IF( TEMP.LE.ZERO )
363 $ ETA = ETA / TWO
364 *
365 TAU = TAU + ETA
366 ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
367 DO 70 J = 1, N
368 DELTA( J ) = DELTA( J ) - ETA
369 WORK( J ) = WORK( J ) + ETA
370 70 CONTINUE
371 *
372 SIGMA = SIGMA + ETA
373 *
374 * Evaluate PSI and the derivative DPSI
375 *
376 DPSI = ZERO
377 PSI = ZERO
378 ERRETM = ZERO
379 DO 80 J = 1, II
380 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
381 PSI = PSI + Z( J )*TEMP
382 DPSI = DPSI + TEMP*TEMP
383 ERRETM = ERRETM + PSI
384 80 CONTINUE
385 ERRETM = ABS( ERRETM )
386 *
387 * Evaluate PHI and the derivative DPHI
388 *
389 TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )
390 PHI = Z( N )*TEMP
391 DPHI = TEMP*TEMP
392 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
393 $ ABS( TAU )*( DPSI+DPHI )
394 *
395 W = RHOINV + PHI + PSI
396 90 CONTINUE
397 *
398 * Return with INFO = 1, NITER = MAXIT and not converged
399 *
400 INFO = 1
401 GO TO 240
402 *
403 * End for the case I = N
404 *
405 ELSE
406 *
407 * The case for I < N
408 *
409 NITER = 1
410 IP1 = I + 1
411 *
412 * Calculate initial guess
413 *
414 DELSQ = ( D( IP1 )-D( I ) )*( D( IP1 )+D( I ) )
415 DELSQ2 = DELSQ / TWO
416 TEMP = DELSQ2 / ( D( I )+SQRT( D( I )*D( I )+DELSQ2 ) )
417 DO 100 J = 1, N
418 WORK( J ) = D( J ) + D( I ) + TEMP
419 DELTA( J ) = ( D( J )-D( I ) ) - TEMP
420 100 CONTINUE
421 *
422 PSI = ZERO
423 DO 110 J = 1, I - 1
424 PSI = PSI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
425 110 CONTINUE
426 *
427 PHI = ZERO
428 DO 120 J = N, I + 2, -1
429 PHI = PHI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
430 120 CONTINUE
431 C = RHOINV + PSI + PHI
432 W = C + Z( I )*Z( I ) / ( WORK( I )*DELTA( I ) ) +
433 $ Z( IP1 )*Z( IP1 ) / ( WORK( IP1 )*DELTA( IP1 ) )
434 *
435 IF( W.GT.ZERO ) THEN
436 *
437 * d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2
438 *
439 * We choose d(i) as origin.
440 *
441 ORGATI = .TRUE.
442 SG2LB = ZERO
443 SG2UB = DELSQ2
444 A = C*DELSQ + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
445 B = Z( I )*Z( I )*DELSQ
446 IF( A.GT.ZERO ) THEN
447 TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
448 ELSE
449 TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
450 END IF
451 *
452 * TAU now is an estimation of SIGMA^2 - D( I )^2. The
453 * following, however, is the corresponding estimation of
454 * SIGMA - D( I ).
455 *
456 ETA = TAU / ( D( I )+SQRT( D( I )*D( I )+TAU ) )
457 ELSE
458 *
459 * (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2
460 *
461 * We choose d(i+1) as origin.
462 *
463 ORGATI = .FALSE.
464 SG2LB = -DELSQ2
465 SG2UB = ZERO
466 A = C*DELSQ - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
467 B = Z( IP1 )*Z( IP1 )*DELSQ
468 IF( A.LT.ZERO ) THEN
469 TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
470 ELSE
471 TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
472 END IF
473 *
474 * TAU now is an estimation of SIGMA^2 - D( IP1 )^2. The
475 * following, however, is the corresponding estimation of
476 * SIGMA - D( IP1 ).
477 *
478 ETA = TAU / ( D( IP1 )+SQRT( ABS( D( IP1 )*D( IP1 )+
479 $ TAU ) ) )
480 END IF
481 *
482 IF( ORGATI ) THEN
483 II = I
484 SIGMA = D( I ) + ETA
485 DO 130 J = 1, N
486 WORK( J ) = D( J ) + D( I ) + ETA
487 DELTA( J ) = ( D( J )-D( I ) ) - ETA
488 130 CONTINUE
489 ELSE
490 II = I + 1
491 SIGMA = D( IP1 ) + ETA
492 DO 140 J = 1, N
493 WORK( J ) = D( J ) + D( IP1 ) + ETA
494 DELTA( J ) = ( D( J )-D( IP1 ) ) - ETA
495 140 CONTINUE
496 END IF
497 IIM1 = II - 1
498 IIP1 = II + 1
499 *
500 * Evaluate PSI and the derivative DPSI
501 *
502 DPSI = ZERO
503 PSI = ZERO
504 ERRETM = ZERO
505 DO 150 J = 1, IIM1
506 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
507 PSI = PSI + Z( J )*TEMP
508 DPSI = DPSI + TEMP*TEMP
509 ERRETM = ERRETM + PSI
510 150 CONTINUE
511 ERRETM = ABS( ERRETM )
512 *
513 * Evaluate PHI and the derivative DPHI
514 *
515 DPHI = ZERO
516 PHI = ZERO
517 DO 160 J = N, IIP1, -1
518 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
519 PHI = PHI + Z( J )*TEMP
520 DPHI = DPHI + TEMP*TEMP
521 ERRETM = ERRETM + PHI
522 160 CONTINUE
523 *
524 W = RHOINV + PHI + PSI
525 *
526 * W is the value of the secular function with
527 * its ii-th element removed.
528 *
529 SWTCH3 = .FALSE.
530 IF( ORGATI ) THEN
531 IF( W.LT.ZERO )
532 $ SWTCH3 = .TRUE.
533 ELSE
534 IF( W.GT.ZERO )
535 $ SWTCH3 = .TRUE.
536 END IF
537 IF( II.EQ.1 .OR. II.EQ.N )
538 $ SWTCH3 = .FALSE.
539 *
540 TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
541 DW = DPSI + DPHI + TEMP*TEMP
542 TEMP = Z( II )*TEMP
543 W = W + TEMP
544 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
545 $ THREE*ABS( TEMP ) + ABS( TAU )*DW
546 *
547 * Test for convergence
548 *
549 IF( ABS( W ).LE.EPS*ERRETM ) THEN
550 GO TO 240
551 END IF
552 *
553 IF( W.LE.ZERO ) THEN
554 SG2LB = MAX( SG2LB, TAU )
555 ELSE
556 SG2UB = MIN( SG2UB, TAU )
557 END IF
558 *
559 * Calculate the new step
560 *
561 NITER = NITER + 1
562 IF( .NOT.SWTCH3 ) THEN
563 DTIPSQ = WORK( IP1 )*DELTA( IP1 )
564 DTISQ = WORK( I )*DELTA( I )
565 IF( ORGATI ) THEN
566 C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
567 ELSE
568 C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
569 END IF
570 A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
571 B = DTIPSQ*DTISQ*W
572 IF( C.EQ.ZERO ) THEN
573 IF( A.EQ.ZERO ) THEN
574 IF( ORGATI ) THEN
575 A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*( DPSI+DPHI )
576 ELSE
577 A = Z( IP1 )*Z( IP1 ) + DTISQ*DTISQ*( DPSI+DPHI )
578 END IF
579 END IF
580 ETA = B / A
581 ELSE IF( A.LE.ZERO ) THEN
582 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
583 ELSE
584 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
585 END IF
586 ELSE
587 *
588 * Interpolation using THREE most relevant poles
589 *
590 DTIIM = WORK( IIM1 )*DELTA( IIM1 )
591 DTIIP = WORK( IIP1 )*DELTA( IIP1 )
592 TEMP = RHOINV + PSI + PHI
593 IF( ORGATI ) THEN
594 TEMP1 = Z( IIM1 ) / DTIIM
595 TEMP1 = TEMP1*TEMP1
596 C = ( TEMP - DTIIP*( DPSI+DPHI ) ) -
597 $ ( D( IIM1 )-D( IIP1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
598 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
599 IF( DPSI.LT.TEMP1 ) THEN
600 ZZ( 3 ) = DTIIP*DTIIP*DPHI
601 ELSE
602 ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
603 END IF
604 ELSE
605 TEMP1 = Z( IIP1 ) / DTIIP
606 TEMP1 = TEMP1*TEMP1
607 C = ( TEMP - DTIIM*( DPSI+DPHI ) ) -
608 $ ( D( IIP1 )-D( IIM1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
609 IF( DPHI.LT.TEMP1 ) THEN
610 ZZ( 1 ) = DTIIM*DTIIM*DPSI
611 ELSE
612 ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
613 END IF
614 ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
615 END IF
616 ZZ( 2 ) = Z( II )*Z( II )
617 DD( 1 ) = DTIIM
618 DD( 2 ) = DELTA( II )*WORK( II )
619 DD( 3 ) = DTIIP
620 CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
621 IF( INFO.NE.0 )
622 $ GO TO 240
623 END IF
624 *
625 * Note, eta should be positive if w is negative, and
626 * eta should be negative otherwise. However,
627 * if for some reason caused by roundoff, eta*w > 0,
628 * we simply use one Newton step instead. This way
629 * will guarantee eta*w < 0.
630 *
631 IF( W*ETA.GE.ZERO )
632 $ ETA = -W / DW
633 IF( ORGATI ) THEN
634 TEMP1 = WORK( I )*DELTA( I )
635 TEMP = ETA - TEMP1
636 ELSE
637 TEMP1 = WORK( IP1 )*DELTA( IP1 )
638 TEMP = ETA - TEMP1
639 END IF
640 IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN
641 IF( W.LT.ZERO ) THEN
642 ETA = ( SG2UB-TAU ) / TWO
643 ELSE
644 ETA = ( SG2LB-TAU ) / TWO
645 END IF
646 END IF
647 *
648 TAU = TAU + ETA
649 ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
650 *
651 PREW = W
652 *
653 SIGMA = SIGMA + ETA
654 DO 170 J = 1, N
655 WORK( J ) = WORK( J ) + ETA
656 DELTA( J ) = DELTA( J ) - ETA
657 170 CONTINUE
658 *
659 * Evaluate PSI and the derivative DPSI
660 *
661 DPSI = ZERO
662 PSI = ZERO
663 ERRETM = ZERO
664 DO 180 J = 1, IIM1
665 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
666 PSI = PSI + Z( J )*TEMP
667 DPSI = DPSI + TEMP*TEMP
668 ERRETM = ERRETM + PSI
669 180 CONTINUE
670 ERRETM = ABS( ERRETM )
671 *
672 * Evaluate PHI and the derivative DPHI
673 *
674 DPHI = ZERO
675 PHI = ZERO
676 DO 190 J = N, IIP1, -1
677 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
678 PHI = PHI + Z( J )*TEMP
679 DPHI = DPHI + TEMP*TEMP
680 ERRETM = ERRETM + PHI
681 190 CONTINUE
682 *
683 TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
684 DW = DPSI + DPHI + TEMP*TEMP
685 TEMP = Z( II )*TEMP
686 W = RHOINV + PHI + PSI + TEMP
687 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
688 $ THREE*ABS( TEMP ) + ABS( TAU )*DW
689 *
690 IF( W.LE.ZERO ) THEN
691 SG2LB = MAX( SG2LB, TAU )
692 ELSE
693 SG2UB = MIN( SG2UB, TAU )
694 END IF
695 *
696 SWTCH = .FALSE.
697 IF( ORGATI ) THEN
698 IF( -W.GT.ABS( PREW ) / TEN )
699 $ SWTCH = .TRUE.
700 ELSE
701 IF( W.GT.ABS( PREW ) / TEN )
702 $ SWTCH = .TRUE.
703 END IF
704 *
705 * Main loop to update the values of the array DELTA and WORK
706 *
707 ITER = NITER + 1
708 *
709 DO 230 NITER = ITER, MAXIT
710 *
711 * Test for convergence
712 *
713 IF( ABS( W ).LE.EPS*ERRETM ) THEN
714 GO TO 240
715 END IF
716 *
717 * Calculate the new step
718 *
719 IF( .NOT.SWTCH3 ) THEN
720 DTIPSQ = WORK( IP1 )*DELTA( IP1 )
721 DTISQ = WORK( I )*DELTA( I )
722 IF( .NOT.SWTCH ) THEN
723 IF( ORGATI ) THEN
724 C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
725 ELSE
726 C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
727 END IF
728 ELSE
729 TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
730 IF( ORGATI ) THEN
731 DPSI = DPSI + TEMP*TEMP
732 ELSE
733 DPHI = DPHI + TEMP*TEMP
734 END IF
735 C = W - DTISQ*DPSI - DTIPSQ*DPHI
736 END IF
737 A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
738 B = DTIPSQ*DTISQ*W
739 IF( C.EQ.ZERO ) THEN
740 IF( A.EQ.ZERO ) THEN
741 IF( .NOT.SWTCH ) THEN
742 IF( ORGATI ) THEN
743 A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*
744 $ ( DPSI+DPHI )
745 ELSE
746 A = Z( IP1 )*Z( IP1 ) +
747 $ DTISQ*DTISQ*( DPSI+DPHI )
748 END IF
749 ELSE
750 A = DTISQ*DTISQ*DPSI + DTIPSQ*DTIPSQ*DPHI
751 END IF
752 END IF
753 ETA = B / A
754 ELSE IF( A.LE.ZERO ) THEN
755 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
756 ELSE
757 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
758 END IF
759 ELSE
760 *
761 * Interpolation using THREE most relevant poles
762 *
763 DTIIM = WORK( IIM1 )*DELTA( IIM1 )
764 DTIIP = WORK( IIP1 )*DELTA( IIP1 )
765 TEMP = RHOINV + PSI + PHI
766 IF( SWTCH ) THEN
767 C = TEMP - DTIIM*DPSI - DTIIP*DPHI
768 ZZ( 1 ) = DTIIM*DTIIM*DPSI
769 ZZ( 3 ) = DTIIP*DTIIP*DPHI
770 ELSE
771 IF( ORGATI ) THEN
772 TEMP1 = Z( IIM1 ) / DTIIM
773 TEMP1 = TEMP1*TEMP1
774 TEMP2 = ( D( IIM1 )-D( IIP1 ) )*
775 $ ( D( IIM1 )+D( IIP1 ) )*TEMP1
776 C = TEMP - DTIIP*( DPSI+DPHI ) - TEMP2
777 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
778 IF( DPSI.LT.TEMP1 ) THEN
779 ZZ( 3 ) = DTIIP*DTIIP*DPHI
780 ELSE
781 ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
782 END IF
783 ELSE
784 TEMP1 = Z( IIP1 ) / DTIIP
785 TEMP1 = TEMP1*TEMP1
786 TEMP2 = ( D( IIP1 )-D( IIM1 ) )*
787 $ ( D( IIM1 )+D( IIP1 ) )*TEMP1
788 C = TEMP - DTIIM*( DPSI+DPHI ) - TEMP2
789 IF( DPHI.LT.TEMP1 ) THEN
790 ZZ( 1 ) = DTIIM*DTIIM*DPSI
791 ELSE
792 ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
793 END IF
794 ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
795 END IF
796 END IF
797 DD( 1 ) = DTIIM
798 DD( 2 ) = DELTA( II )*WORK( II )
799 DD( 3 ) = DTIIP
800 CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
801 IF( INFO.NE.0 )
802 $ GO TO 240
803 END IF
804 *
805 * Note, eta should be positive if w is negative, and
806 * eta should be negative otherwise. However,
807 * if for some reason caused by roundoff, eta*w > 0,
808 * we simply use one Newton step instead. This way
809 * will guarantee eta*w < 0.
810 *
811 IF( W*ETA.GE.ZERO )
812 $ ETA = -W / DW
813 IF( ORGATI ) THEN
814 TEMP1 = WORK( I )*DELTA( I )
815 TEMP = ETA - TEMP1
816 ELSE
817 TEMP1 = WORK( IP1 )*DELTA( IP1 )
818 TEMP = ETA - TEMP1
819 END IF
820 IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN
821 IF( W.LT.ZERO ) THEN
822 ETA = ( SG2UB-TAU ) / TWO
823 ELSE
824 ETA = ( SG2LB-TAU ) / TWO
825 END IF
826 END IF
827 *
828 TAU = TAU + ETA
829 ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
830 *
831 SIGMA = SIGMA + ETA
832 DO 200 J = 1, N
833 WORK( J ) = WORK( J ) + ETA
834 DELTA( J ) = DELTA( J ) - ETA
835 200 CONTINUE
836 *
837 PREW = W
838 *
839 * Evaluate PSI and the derivative DPSI
840 *
841 DPSI = ZERO
842 PSI = ZERO
843 ERRETM = ZERO
844 DO 210 J = 1, IIM1
845 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
846 PSI = PSI + Z( J )*TEMP
847 DPSI = DPSI + TEMP*TEMP
848 ERRETM = ERRETM + PSI
849 210 CONTINUE
850 ERRETM = ABS( ERRETM )
851 *
852 * Evaluate PHI and the derivative DPHI
853 *
854 DPHI = ZERO
855 PHI = ZERO
856 DO 220 J = N, IIP1, -1
857 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
858 PHI = PHI + Z( J )*TEMP
859 DPHI = DPHI + TEMP*TEMP
860 ERRETM = ERRETM + PHI
861 220 CONTINUE
862 *
863 TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
864 DW = DPSI + DPHI + TEMP*TEMP
865 TEMP = Z( II )*TEMP
866 W = RHOINV + PHI + PSI + TEMP
867 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
868 $ THREE*ABS( TEMP ) + ABS( TAU )*DW
869 IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )
870 $ SWTCH = .NOT.SWTCH
871 *
872 IF( W.LE.ZERO ) THEN
873 SG2LB = MAX( SG2LB, TAU )
874 ELSE
875 SG2UB = MIN( SG2UB, TAU )
876 END IF
877 *
878 230 CONTINUE
879 *
880 * Return with INFO = 1, NITER = MAXIT and not converged
881 *
882 INFO = 1
883 *
884 END IF
885 *
886 240 CONTINUE
887 RETURN
888 *
889 * End of DLASD4
890 *
891 END