1: SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
2: *
3: * -- LAPACK auxiliary routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
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 = 20 )
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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd4.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack