Annotation of rpl/lapack/lapack/dlasd4.f, revision 1.21
1.13 bertrand 1: *> \brief \b DLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. Used by dbdsdc.
1.10 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.13 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.10 bertrand 7: *
8: *> \htmlonly
1.13 bertrand 9: *> Download DLASD4 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd4.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd4.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd4.f">
1.10 bertrand 15: *> [TXT]</a>
1.13 bertrand 16: *> \endhtmlonly
1.10 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
1.13 bertrand 22: *
1.10 bertrand 23: * .. Scalar Arguments ..
24: * INTEGER I, INFO, N
25: * DOUBLE PRECISION RHO, SIGMA
26: * ..
27: * .. Array Arguments ..
28: * DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), Z( * )
29: * ..
1.13 bertrand 30: *
1.10 bertrand 31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> This subroutine computes the square root of the I-th updated
38: *> eigenvalue of a positive symmetric rank-one modification to
39: *> a positive diagonal matrix whose entries are given as the squares
40: *> of the corresponding entries in the array d, and that
41: *>
42: *> 0 <= D(i) < D(j) for i < j
43: *>
44: *> and that RHO > 0. This is arranged by the calling routine, and is
45: *> no loss in generality. The rank-one modified system is thus
46: *>
47: *> diag( D ) * diag( D ) + RHO * Z * Z_transpose.
48: *>
49: *> where we assume the Euclidean norm of Z is 1.
50: *>
51: *> The method consists of approximating the rational functions in the
52: *> secular equation by simpler interpolating rational functions.
53: *> \endverbatim
54: *
55: * Arguments:
56: * ==========
57: *
58: *> \param[in] N
59: *> \verbatim
60: *> N is INTEGER
61: *> The length of all arrays.
62: *> \endverbatim
63: *>
64: *> \param[in] I
65: *> \verbatim
66: *> I is INTEGER
67: *> The index of the eigenvalue to be computed. 1 <= I <= N.
68: *> \endverbatim
69: *>
70: *> \param[in] D
71: *> \verbatim
72: *> D is DOUBLE PRECISION array, dimension ( N )
73: *> The original eigenvalues. It is assumed that they are in
74: *> order, 0 <= D(I) < D(J) for I < J.
75: *> \endverbatim
76: *>
77: *> \param[in] Z
78: *> \verbatim
79: *> Z is DOUBLE PRECISION array, dimension ( N )
80: *> The components of the updating vector.
81: *> \endverbatim
82: *>
83: *> \param[out] DELTA
84: *> \verbatim
85: *> DELTA is DOUBLE PRECISION array, dimension ( N )
86: *> If N .ne. 1, DELTA contains (D(j) - sigma_I) in its j-th
87: *> component. If N = 1, then DELTA(1) = 1. The vector DELTA
88: *> contains the information necessary to construct the
89: *> (singular) eigenvectors.
90: *> \endverbatim
91: *>
92: *> \param[in] RHO
93: *> \verbatim
94: *> RHO is DOUBLE PRECISION
95: *> The scalar in the symmetric updating formula.
96: *> \endverbatim
97: *>
98: *> \param[out] SIGMA
99: *> \verbatim
100: *> SIGMA is DOUBLE PRECISION
101: *> The computed sigma_I, the I-th updated eigenvalue.
102: *> \endverbatim
103: *>
104: *> \param[out] WORK
105: *> \verbatim
106: *> WORK is DOUBLE PRECISION array, dimension ( N )
107: *> If N .ne. 1, WORK contains (D(j) + sigma_I) in its j-th
108: *> component. If N = 1, then WORK( 1 ) = 1.
109: *> \endverbatim
110: *>
111: *> \param[out] INFO
112: *> \verbatim
113: *> INFO is INTEGER
114: *> = 0: successful exit
115: *> > 0: if INFO = 1, the updating process failed.
116: *> \endverbatim
117: *
118: *> \par Internal Parameters:
119: * =========================
120: *>
121: *> \verbatim
122: *> Logical variable ORGATI (origin-at-i?) is used for distinguishing
123: *> whether D(i) or D(i+1) is treated as the origin.
124: *>
125: *> ORGATI = .true. origin at i
126: *> ORGATI = .false. origin at i+1
127: *>
128: *> Logical variable SWTCH3 (switch-for-3-poles?) is for noting
129: *> if we are working with THREE poles!
130: *>
131: *> MAXIT is the maximum number of iterations allowed for each
132: *> eigenvalue.
133: *> \endverbatim
134: *
135: * Authors:
136: * ========
137: *
1.13 bertrand 138: *> \author Univ. of Tennessee
139: *> \author Univ. of California Berkeley
140: *> \author Univ. of Colorado Denver
141: *> \author NAG Ltd.
1.10 bertrand 142: *
1.18 bertrand 143: *> \ingroup OTHERauxiliary
1.10 bertrand 144: *
145: *> \par Contributors:
146: * ==================
147: *>
148: *> Ren-Cang Li, Computer Science Division, University of California
149: *> at Berkeley, USA
150: *>
151: * =====================================================================
1.1 bertrand 152: SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
153: *
1.21 ! bertrand 154: * -- LAPACK auxiliary routine --
1.1 bertrand 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: DOUBLE PRECISION RHO, SIGMA
161: * ..
162: * .. Array Arguments ..
163: DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), Z( * )
164: * ..
165: *
166: * =====================================================================
167: *
168: * .. Parameters ..
169: INTEGER MAXIT
1.13 bertrand 170: PARAMETER ( MAXIT = 400 )
1.1 bertrand 171: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
172: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
173: $ THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0,
174: $ TEN = 10.0D+0 )
175: * ..
176: * .. Local Scalars ..
1.13 bertrand 177: LOGICAL ORGATI, SWTCH, SWTCH3, GEOMAVG
1.1 bertrand 178: INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
1.13 bertrand 179: DOUBLE PRECISION A, B, C, DELSQ, DELSQ2, SQ2, DPHI, DPSI, DTIIM,
1.1 bertrand 180: $ DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS,
1.13 bertrand 181: $ ERRETM, ETA, PHI, PREW, PSI, RHOINV, SGLB,
182: $ SGUB, TAU, TAU2, TEMP, TEMP1, TEMP2, W
1.1 bertrand 183: * ..
184: * .. Local Arrays ..
185: DOUBLE PRECISION DD( 3 ), ZZ( 3 )
186: * ..
187: * .. External Subroutines ..
188: EXTERNAL DLAED6, DLASD5
189: * ..
190: * .. External Functions ..
191: DOUBLE PRECISION DLAMCH
192: EXTERNAL DLAMCH
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 DLASD5( I, D, Z, DELTA, RHO, SIGMA, WORK )
216: RETURN
217: END IF
218: *
219: * Compute machine epsilon
220: *
221: EPS = DLAMCH( 'Epsilon' )
222: RHOINV = ONE / RHO
1.15 bertrand 223: TAU2= ZERO
1.1 bertrand 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: *
1.13 bertrand 262: * The following TAU2 is to approximate
1.1 bertrand 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
1.13 bertrand 272: TAU2 = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
1.1 bertrand 273: ELSE
1.13 bertrand 274: TAU2 = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
1.1 bertrand 275: END IF
1.15 bertrand 276: TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
1.1 bertrand 277: END IF
278: *
279: * It can be proved that
1.13 bertrand 280: * D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU2 <= D(N)^2+RHO
1.1 bertrand 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: *
1.13 bertrand 287: * The following TAU2 is to approximate
1.1 bertrand 288: * SIGMA_n^2 - D( N )*D( N )
289: *
290: IF( A.LT.ZERO ) THEN
1.13 bertrand 291: TAU2 = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
1.1 bertrand 292: ELSE
1.13 bertrand 293: TAU2 = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
1.1 bertrand 294: END IF
1.15 bertrand 295: TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
296:
1.1 bertrand 297: *
298: * It can be proved that
1.13 bertrand 299: * D(N)^2 < D(N)^2+TAU2 < SIGMA(N)^2 < D(N)^2+RHO/2
1.1 bertrand 300: *
301: END IF
302: *
1.13 bertrand 303: * The following TAU is to approximate SIGMA_n - D( N )
1.1 bertrand 304: *
1.15 bertrand 305: * TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
1.1 bertrand 306: *
1.13 bertrand 307: SIGMA = D( N ) + TAU
1.1 bertrand 308: DO 30 J = 1, N
1.13 bertrand 309: DELTA( J ) = ( D( J )-D( N ) ) - TAU
310: WORK( J ) = D( J ) + D( N ) + TAU
1.1 bertrand 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
1.18 bertrand 331: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
1.13 bertrand 332: * $ + ABS( TAU2 )*( DPSI+DPHI )
1.1 bertrand 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: *
1.13 bertrand 372: ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
1.1 bertrand 373: TAU = TAU + ETA
1.13 bertrand 374: SIGMA = SIGMA + ETA
375: *
1.1 bertrand 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: *
1.13 bertrand 396: TAU2 = WORK( N )*DELTA( N )
397: TEMP = Z( N ) / TAU2
1.1 bertrand 398: PHI = Z( N )*TEMP
399: DPHI = TEMP*TEMP
1.18 bertrand 400: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
1.13 bertrand 401: * $ + ABS( TAU2 )*( DPSI+DPHI )
1.1 bertrand 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: *
1.13 bertrand 442: ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
1.1 bertrand 443: TAU = TAU + ETA
1.13 bertrand 444: SIGMA = SIGMA + ETA
445: *
1.1 bertrand 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: *
1.13 bertrand 466: TAU2 = WORK( N )*DELTA( N )
467: TEMP = Z( N ) / TAU2
1.1 bertrand 468: PHI = Z( N )*TEMP
469: DPHI = TEMP*TEMP
1.18 bertrand 470: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
1.13 bertrand 471: * $ + ABS( TAU2 )*( DPSI+DPHI )
1.1 bertrand 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
1.13 bertrand 494: SQ2=SQRT( ( D( I )*D( I )+D( IP1 )*D( IP1 ) ) / TWO )
495: TEMP = DELSQ2 / ( D( I )+SQ2 )
1.1 bertrand 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: *
1.13 bertrand 514: GEOMAVG = .FALSE.
1.1 bertrand 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.
1.13 bertrand 522: II = I
523: SGLB = ZERO
524: SGUB = DELSQ2 / ( D( I )+SQ2 )
1.1 bertrand 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
1.13 bertrand 528: TAU2 = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
1.1 bertrand 529: ELSE
1.13 bertrand 530: TAU2 = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
1.1 bertrand 531: END IF
532: *
1.13 bertrand 533: * TAU2 now is an estimation of SIGMA^2 - D( I )^2. The
1.1 bertrand 534: * following, however, is the corresponding estimation of
535: * SIGMA - D( I ).
536: *
1.13 bertrand 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
1.1 bertrand 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.
1.13 bertrand 551: II = IP1
552: SGLB = -DELSQ2 / ( D( II )+SQ2 )
553: SGUB = ZERO
1.1 bertrand 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
1.13 bertrand 557: TAU2 = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
1.1 bertrand 558: ELSE
1.13 bertrand 559: TAU2 = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
1.1 bertrand 560: END IF
561: *
1.13 bertrand 562: * TAU2 now is an estimation of SIGMA^2 - D( IP1 )^2. The
1.1 bertrand 563: * following, however, is the corresponding estimation of
564: * SIGMA - D( IP1 ).
565: *
1.13 bertrand 566: TAU = TAU2 / ( D( IP1 )+SQRT( ABS( D( IP1 )*D( IP1 )+
567: $ TAU2 ) ) )
1.1 bertrand 568: END IF
569: *
1.13 bertrand 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
1.1 bertrand 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
1.18 bertrand 622: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
623: $ + THREE*ABS( TEMP )
1.13 bertrand 624: * $ + ABS( TAU2 )*DW
1.1 bertrand 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
1.13 bertrand 633: SGLB = MAX( SGLB, TAU )
1.1 bertrand 634: ELSE
1.13 bertrand 635: SGUB = MIN( SGUB, TAU )
1.1 bertrand 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 DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
1.13 bertrand 700: *
701: IF( INFO.NE.0 ) THEN
702: *
1.18 bertrand 703: * If INFO is not 0, i.e., DLAED6 failed, switch back
1.13 bertrand 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
1.1 bertrand 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
1.13 bertrand 742: *
743: ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
744: TEMP = TAU + ETA
745: IF( TEMP.GT.SGUB .OR. TEMP.LT.SGLB ) THEN
1.1 bertrand 746: IF( W.LT.ZERO ) THEN
1.13 bertrand 747: ETA = ( SGUB-TAU ) / TWO
1.1 bertrand 748: ELSE
1.13 bertrand 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
1.1 bertrand 761: END IF
762: END IF
763: *
764: PREW = W
765: *
1.13 bertrand 766: TAU = TAU + ETA
1.1 bertrand 767: SIGMA = SIGMA + ETA
1.13 bertrand 768: *
1.1 bertrand 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: *
1.13 bertrand 798: TAU2 = WORK( II )*DELTA( II )
799: TEMP = Z( II ) / TAU2
1.1 bertrand 800: DW = DPSI + DPHI + TEMP*TEMP
801: TEMP = Z( II )*TEMP
802: W = RHOINV + PHI + PSI + TEMP
1.18 bertrand 803: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
804: $ + THREE*ABS( TEMP )
1.13 bertrand 805: * $ + ABS( TAU2 )*DW
1.1 bertrand 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
1.13 bertrand 825: * $ .OR. (SGUB-SGLB).LE.EIGHT*ABS(SGUB+SGLB) ) THEN
1.1 bertrand 826: GO TO 240
827: END IF
828: *
1.13 bertrand 829: IF( W.LE.ZERO ) THEN
830: SGLB = MAX( SGLB, TAU )
831: ELSE
832: SGUB = MIN( SGUB, TAU )
833: END IF
834: *
1.1 bertrand 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 DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
1.13 bertrand 919: *
920: IF( INFO.NE.0 ) THEN
921: *
1.18 bertrand 922: * If INFO is not 0, i.e., DLAED6 failed, switch
1.13 bertrand 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
1.1 bertrand 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
1.13 bertrand 977: *
978: ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
979: TEMP=TAU+ETA
980: IF( TEMP.GT.SGUB .OR. TEMP.LT.SGLB ) THEN
1.1 bertrand 981: IF( W.LT.ZERO ) THEN
1.13 bertrand 982: ETA = ( SGUB-TAU ) / TWO
1.1 bertrand 983: ELSE
1.13 bertrand 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
1.1 bertrand 996: END IF
997: END IF
998: *
1.13 bertrand 999: PREW = W
1000: *
1.1 bertrand 1001: TAU = TAU + ETA
1.13 bertrand 1002: SIGMA = SIGMA + ETA
1.1 bertrand 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: *
1.13 bertrand 1033: TAU2 = WORK( II )*DELTA( II )
1034: TEMP = Z( II ) / TAU2
1.1 bertrand 1035: DW = DPSI + DPHI + TEMP*TEMP
1036: TEMP = Z( II )*TEMP
1037: W = RHOINV + PHI + PSI + TEMP
1.18 bertrand 1038: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
1039: $ + THREE*ABS( TEMP )
1.13 bertrand 1040: * $ + ABS( TAU2 )*DW
1041: *
1.1 bertrand 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 DLASD4
1057: *
1058: 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