Annotation of rpl/lapack/lapack/dlasd4.f, revision 1.17
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.15 bertrand 143: *> \date November 2013
1.10 bertrand 144: *
145: *> \ingroup auxOTHERauxiliary
146: *
147: *> \par Contributors:
148: * ==================
149: *>
150: *> Ren-Cang Li, Computer Science Division, University of California
151: *> at Berkeley, USA
152: *>
153: * =====================================================================
1.1 bertrand 154: SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
155: *
1.15 bertrand 156: * -- LAPACK auxiliary routine (version 3.5.0) --
1.1 bertrand 157: * -- LAPACK is a software package provided by Univ. of Tennessee, --
158: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.15 bertrand 159: * November 2013
1.1 bertrand 160: *
161: * .. Scalar Arguments ..
162: INTEGER I, INFO, N
163: DOUBLE PRECISION RHO, SIGMA
164: * ..
165: * .. Array Arguments ..
166: DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), Z( * )
167: * ..
168: *
169: * =====================================================================
170: *
171: * .. Parameters ..
172: INTEGER MAXIT
1.13 bertrand 173: PARAMETER ( MAXIT = 400 )
1.1 bertrand 174: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
175: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
176: $ THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0,
177: $ TEN = 10.0D+0 )
178: * ..
179: * .. Local Scalars ..
1.13 bertrand 180: LOGICAL ORGATI, SWTCH, SWTCH3, GEOMAVG
1.1 bertrand 181: INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
1.13 bertrand 182: DOUBLE PRECISION A, B, C, DELSQ, DELSQ2, SQ2, DPHI, DPSI, DTIIM,
1.1 bertrand 183: $ DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS,
1.13 bertrand 184: $ ERRETM, ETA, PHI, PREW, PSI, RHOINV, SGLB,
185: $ SGUB, TAU, TAU2, TEMP, TEMP1, TEMP2, W
1.1 bertrand 186: * ..
187: * .. Local Arrays ..
188: DOUBLE PRECISION DD( 3 ), ZZ( 3 )
189: * ..
190: * .. External Subroutines ..
191: EXTERNAL DLAED6, DLASD5
192: * ..
193: * .. External Functions ..
194: DOUBLE PRECISION DLAMCH
195: EXTERNAL DLAMCH
196: * ..
197: * .. Intrinsic Functions ..
198: INTRINSIC ABS, MAX, MIN, SQRT
199: * ..
200: * .. Executable Statements ..
201: *
202: * Since this routine is called in an inner loop, we do no argument
203: * checking.
204: *
205: * Quick return for N=1 and 2.
206: *
207: INFO = 0
208: IF( N.EQ.1 ) THEN
209: *
210: * Presumably, I=1 upon entry
211: *
212: SIGMA = SQRT( D( 1 )*D( 1 )+RHO*Z( 1 )*Z( 1 ) )
213: DELTA( 1 ) = ONE
214: WORK( 1 ) = ONE
215: RETURN
216: END IF
217: IF( N.EQ.2 ) THEN
218: CALL DLASD5( I, D, Z, DELTA, RHO, SIGMA, WORK )
219: RETURN
220: END IF
221: *
222: * Compute machine epsilon
223: *
224: EPS = DLAMCH( 'Epsilon' )
225: RHOINV = ONE / RHO
1.15 bertrand 226: TAU2= ZERO
1.1 bertrand 227: *
228: * The case I = N
229: *
230: IF( I.EQ.N ) THEN
231: *
232: * Initialize some basic variables
233: *
234: II = N - 1
235: NITER = 1
236: *
237: * Calculate initial guess
238: *
239: TEMP = RHO / TWO
240: *
241: * If ||Z||_2 is not one, then TEMP should be set to
242: * RHO * ||Z||_2^2 / TWO
243: *
244: TEMP1 = TEMP / ( D( N )+SQRT( D( N )*D( N )+TEMP ) )
245: DO 10 J = 1, N
246: WORK( J ) = D( J ) + D( N ) + TEMP1
247: DELTA( J ) = ( D( J )-D( N ) ) - TEMP1
248: 10 CONTINUE
249: *
250: PSI = ZERO
251: DO 20 J = 1, N - 2
252: PSI = PSI + Z( J )*Z( J ) / ( DELTA( J )*WORK( J ) )
253: 20 CONTINUE
254: *
255: C = RHOINV + PSI
256: W = C + Z( II )*Z( II ) / ( DELTA( II )*WORK( II ) ) +
257: $ Z( N )*Z( N ) / ( DELTA( N )*WORK( N ) )
258: *
259: IF( W.LE.ZERO ) THEN
260: TEMP1 = SQRT( D( N )*D( N )+RHO )
261: TEMP = Z( N-1 )*Z( N-1 ) / ( ( D( N-1 )+TEMP1 )*
262: $ ( D( N )-D( N-1 )+RHO / ( D( N )+TEMP1 ) ) ) +
263: $ Z( N )*Z( N ) / RHO
264: *
1.13 bertrand 265: * The following TAU2 is to approximate
1.1 bertrand 266: * SIGMA_n^2 - D( N )*D( N )
267: *
268: IF( C.LE.TEMP ) THEN
269: TAU = RHO
270: ELSE
271: DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
272: A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
273: B = Z( N )*Z( N )*DELSQ
274: IF( A.LT.ZERO ) THEN
1.13 bertrand 275: TAU2 = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
1.1 bertrand 276: ELSE
1.13 bertrand 277: TAU2 = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
1.1 bertrand 278: END IF
1.15 bertrand 279: TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
1.1 bertrand 280: END IF
281: *
282: * It can be proved that
1.13 bertrand 283: * D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU2 <= D(N)^2+RHO
1.1 bertrand 284: *
285: ELSE
286: DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
287: A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
288: B = Z( N )*Z( N )*DELSQ
289: *
1.13 bertrand 290: * The following TAU2 is to approximate
1.1 bertrand 291: * SIGMA_n^2 - D( N )*D( N )
292: *
293: IF( A.LT.ZERO ) THEN
1.13 bertrand 294: TAU2 = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
1.1 bertrand 295: ELSE
1.13 bertrand 296: TAU2 = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
1.1 bertrand 297: END IF
1.15 bertrand 298: TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
299:
1.1 bertrand 300: *
301: * It can be proved that
1.13 bertrand 302: * D(N)^2 < D(N)^2+TAU2 < SIGMA(N)^2 < D(N)^2+RHO/2
1.1 bertrand 303: *
304: END IF
305: *
1.13 bertrand 306: * The following TAU is to approximate SIGMA_n - D( N )
1.1 bertrand 307: *
1.15 bertrand 308: * TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
1.1 bertrand 309: *
1.13 bertrand 310: SIGMA = D( N ) + TAU
1.1 bertrand 311: DO 30 J = 1, N
1.13 bertrand 312: DELTA( J ) = ( D( J )-D( N ) ) - TAU
313: WORK( J ) = D( J ) + D( N ) + TAU
1.1 bertrand 314: 30 CONTINUE
315: *
316: * Evaluate PSI and the derivative DPSI
317: *
318: DPSI = ZERO
319: PSI = ZERO
320: ERRETM = ZERO
321: DO 40 J = 1, II
322: TEMP = Z( J ) / ( DELTA( J )*WORK( J ) )
323: PSI = PSI + Z( J )*TEMP
324: DPSI = DPSI + TEMP*TEMP
325: ERRETM = ERRETM + PSI
326: 40 CONTINUE
327: ERRETM = ABS( ERRETM )
328: *
329: * Evaluate PHI and the derivative DPHI
330: *
331: TEMP = Z( N ) / ( DELTA( N )*WORK( N ) )
332: PHI = Z( N )*TEMP
333: DPHI = TEMP*TEMP
1.13 bertrand 334: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
335: * $ + ABS( TAU2 )*( DPSI+DPHI )
1.1 bertrand 336: *
337: W = RHOINV + PHI + PSI
338: *
339: * Test for convergence
340: *
341: IF( ABS( W ).LE.EPS*ERRETM ) THEN
342: GO TO 240
343: END IF
344: *
345: * Calculate the new step
346: *
347: NITER = NITER + 1
348: DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
349: DTNSQ = WORK( N )*DELTA( N )
350: C = W - DTNSQ1*DPSI - DTNSQ*DPHI
351: A = ( DTNSQ+DTNSQ1 )*W - DTNSQ*DTNSQ1*( DPSI+DPHI )
352: B = DTNSQ*DTNSQ1*W
353: IF( C.LT.ZERO )
354: $ C = ABS( C )
355: IF( C.EQ.ZERO ) THEN
356: ETA = RHO - SIGMA*SIGMA
357: ELSE IF( A.GE.ZERO ) THEN
358: ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
359: ELSE
360: ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
361: END IF
362: *
363: * Note, eta should be positive if w is negative, and
364: * eta should be negative otherwise. However,
365: * if for some reason caused by roundoff, eta*w > 0,
366: * we simply use one Newton step instead. This way
367: * will guarantee eta*w < 0.
368: *
369: IF( W*ETA.GT.ZERO )
370: $ ETA = -W / ( DPSI+DPHI )
371: TEMP = ETA - DTNSQ
372: IF( TEMP.GT.RHO )
373: $ ETA = RHO + DTNSQ
374: *
1.13 bertrand 375: ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
1.1 bertrand 376: TAU = TAU + ETA
1.13 bertrand 377: SIGMA = SIGMA + ETA
378: *
1.1 bertrand 379: DO 50 J = 1, N
380: DELTA( J ) = DELTA( J ) - ETA
381: WORK( J ) = WORK( J ) + ETA
382: 50 CONTINUE
383: *
384: * Evaluate PSI and the derivative DPSI
385: *
386: DPSI = ZERO
387: PSI = ZERO
388: ERRETM = ZERO
389: DO 60 J = 1, II
390: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
391: PSI = PSI + Z( J )*TEMP
392: DPSI = DPSI + TEMP*TEMP
393: ERRETM = ERRETM + PSI
394: 60 CONTINUE
395: ERRETM = ABS( ERRETM )
396: *
397: * Evaluate PHI and the derivative DPHI
398: *
1.13 bertrand 399: TAU2 = WORK( N )*DELTA( N )
400: TEMP = Z( N ) / TAU2
1.1 bertrand 401: PHI = Z( N )*TEMP
402: DPHI = TEMP*TEMP
1.13 bertrand 403: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
404: * $ + ABS( TAU2 )*( DPSI+DPHI )
1.1 bertrand 405: *
406: W = RHOINV + PHI + PSI
407: *
408: * Main loop to update the values of the array DELTA
409: *
410: ITER = NITER + 1
411: *
412: DO 90 NITER = ITER, MAXIT
413: *
414: * Test for convergence
415: *
416: IF( ABS( W ).LE.EPS*ERRETM ) THEN
417: GO TO 240
418: END IF
419: *
420: * Calculate the new step
421: *
422: DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
423: DTNSQ = WORK( N )*DELTA( N )
424: C = W - DTNSQ1*DPSI - DTNSQ*DPHI
425: A = ( DTNSQ+DTNSQ1 )*W - DTNSQ1*DTNSQ*( DPSI+DPHI )
426: B = DTNSQ1*DTNSQ*W
427: IF( A.GE.ZERO ) THEN
428: ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
429: ELSE
430: ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
431: END IF
432: *
433: * Note, eta should be positive if w is negative, and
434: * eta should be negative otherwise. However,
435: * if for some reason caused by roundoff, eta*w > 0,
436: * we simply use one Newton step instead. This way
437: * will guarantee eta*w < 0.
438: *
439: IF( W*ETA.GT.ZERO )
440: $ ETA = -W / ( DPSI+DPHI )
441: TEMP = ETA - DTNSQ
442: IF( TEMP.LE.ZERO )
443: $ ETA = ETA / TWO
444: *
1.13 bertrand 445: ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
1.1 bertrand 446: TAU = TAU + ETA
1.13 bertrand 447: SIGMA = SIGMA + ETA
448: *
1.1 bertrand 449: DO 70 J = 1, N
450: DELTA( J ) = DELTA( J ) - ETA
451: WORK( J ) = WORK( J ) + ETA
452: 70 CONTINUE
453: *
454: * Evaluate PSI and the derivative DPSI
455: *
456: DPSI = ZERO
457: PSI = ZERO
458: ERRETM = ZERO
459: DO 80 J = 1, II
460: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
461: PSI = PSI + Z( J )*TEMP
462: DPSI = DPSI + TEMP*TEMP
463: ERRETM = ERRETM + PSI
464: 80 CONTINUE
465: ERRETM = ABS( ERRETM )
466: *
467: * Evaluate PHI and the derivative DPHI
468: *
1.13 bertrand 469: TAU2 = WORK( N )*DELTA( N )
470: TEMP = Z( N ) / TAU2
1.1 bertrand 471: PHI = Z( N )*TEMP
472: DPHI = TEMP*TEMP
1.13 bertrand 473: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
474: * $ + ABS( TAU2 )*( DPSI+DPHI )
1.1 bertrand 475: *
476: W = RHOINV + PHI + PSI
477: 90 CONTINUE
478: *
479: * Return with INFO = 1, NITER = MAXIT and not converged
480: *
481: INFO = 1
482: GO TO 240
483: *
484: * End for the case I = N
485: *
486: ELSE
487: *
488: * The case for I < N
489: *
490: NITER = 1
491: IP1 = I + 1
492: *
493: * Calculate initial guess
494: *
495: DELSQ = ( D( IP1 )-D( I ) )*( D( IP1 )+D( I ) )
496: DELSQ2 = DELSQ / TWO
1.13 bertrand 497: SQ2=SQRT( ( D( I )*D( I )+D( IP1 )*D( IP1 ) ) / TWO )
498: TEMP = DELSQ2 / ( D( I )+SQ2 )
1.1 bertrand 499: DO 100 J = 1, N
500: WORK( J ) = D( J ) + D( I ) + TEMP
501: DELTA( J ) = ( D( J )-D( I ) ) - TEMP
502: 100 CONTINUE
503: *
504: PSI = ZERO
505: DO 110 J = 1, I - 1
506: PSI = PSI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
507: 110 CONTINUE
508: *
509: PHI = ZERO
510: DO 120 J = N, I + 2, -1
511: PHI = PHI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
512: 120 CONTINUE
513: C = RHOINV + PSI + PHI
514: W = C + Z( I )*Z( I ) / ( WORK( I )*DELTA( I ) ) +
515: $ Z( IP1 )*Z( IP1 ) / ( WORK( IP1 )*DELTA( IP1 ) )
516: *
1.13 bertrand 517: GEOMAVG = .FALSE.
1.1 bertrand 518: IF( W.GT.ZERO ) THEN
519: *
520: * d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2
521: *
522: * We choose d(i) as origin.
523: *
524: ORGATI = .TRUE.
1.13 bertrand 525: II = I
526: SGLB = ZERO
527: SGUB = DELSQ2 / ( D( I )+SQ2 )
1.1 bertrand 528: A = C*DELSQ + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
529: B = Z( I )*Z( I )*DELSQ
530: IF( A.GT.ZERO ) THEN
1.13 bertrand 531: TAU2 = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
1.1 bertrand 532: ELSE
1.13 bertrand 533: TAU2 = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
1.1 bertrand 534: END IF
535: *
1.13 bertrand 536: * TAU2 now is an estimation of SIGMA^2 - D( I )^2. The
1.1 bertrand 537: * following, however, is the corresponding estimation of
538: * SIGMA - D( I ).
539: *
1.13 bertrand 540: TAU = TAU2 / ( D( I )+SQRT( D( I )*D( I )+TAU2 ) )
541: TEMP = SQRT(EPS)
542: IF( (D(I).LE.TEMP*D(IP1)).AND.(ABS(Z(I)).LE.TEMP)
543: $ .AND.(D(I).GT.ZERO) ) THEN
544: TAU = MIN( TEN*D(I), SGUB )
545: GEOMAVG = .TRUE.
546: END IF
1.1 bertrand 547: ELSE
548: *
549: * (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2
550: *
551: * We choose d(i+1) as origin.
552: *
553: ORGATI = .FALSE.
1.13 bertrand 554: II = IP1
555: SGLB = -DELSQ2 / ( D( II )+SQ2 )
556: SGUB = ZERO
1.1 bertrand 557: A = C*DELSQ - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
558: B = Z( IP1 )*Z( IP1 )*DELSQ
559: IF( A.LT.ZERO ) THEN
1.13 bertrand 560: TAU2 = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
1.1 bertrand 561: ELSE
1.13 bertrand 562: TAU2 = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
1.1 bertrand 563: END IF
564: *
1.13 bertrand 565: * TAU2 now is an estimation of SIGMA^2 - D( IP1 )^2. The
1.1 bertrand 566: * following, however, is the corresponding estimation of
567: * SIGMA - D( IP1 ).
568: *
1.13 bertrand 569: TAU = TAU2 / ( D( IP1 )+SQRT( ABS( D( IP1 )*D( IP1 )+
570: $ TAU2 ) ) )
1.1 bertrand 571: END IF
572: *
1.13 bertrand 573: SIGMA = D( II ) + TAU
574: DO 130 J = 1, N
575: WORK( J ) = D( J ) + D( II ) + TAU
576: DELTA( J ) = ( D( J )-D( II ) ) - TAU
577: 130 CONTINUE
1.1 bertrand 578: IIM1 = II - 1
579: IIP1 = II + 1
580: *
581: * Evaluate PSI and the derivative DPSI
582: *
583: DPSI = ZERO
584: PSI = ZERO
585: ERRETM = ZERO
586: DO 150 J = 1, IIM1
587: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
588: PSI = PSI + Z( J )*TEMP
589: DPSI = DPSI + TEMP*TEMP
590: ERRETM = ERRETM + PSI
591: 150 CONTINUE
592: ERRETM = ABS( ERRETM )
593: *
594: * Evaluate PHI and the derivative DPHI
595: *
596: DPHI = ZERO
597: PHI = ZERO
598: DO 160 J = N, IIP1, -1
599: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
600: PHI = PHI + Z( J )*TEMP
601: DPHI = DPHI + TEMP*TEMP
602: ERRETM = ERRETM + PHI
603: 160 CONTINUE
604: *
605: W = RHOINV + PHI + PSI
606: *
607: * W is the value of the secular function with
608: * its ii-th element removed.
609: *
610: SWTCH3 = .FALSE.
611: IF( ORGATI ) THEN
612: IF( W.LT.ZERO )
613: $ SWTCH3 = .TRUE.
614: ELSE
615: IF( W.GT.ZERO )
616: $ SWTCH3 = .TRUE.
617: END IF
618: IF( II.EQ.1 .OR. II.EQ.N )
619: $ SWTCH3 = .FALSE.
620: *
621: TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
622: DW = DPSI + DPHI + TEMP*TEMP
623: TEMP = Z( II )*TEMP
624: W = W + TEMP
1.13 bertrand 625: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
626: $ + THREE*ABS( TEMP )
627: * $ + ABS( TAU2 )*DW
1.1 bertrand 628: *
629: * Test for convergence
630: *
631: IF( ABS( W ).LE.EPS*ERRETM ) THEN
632: GO TO 240
633: END IF
634: *
635: IF( W.LE.ZERO ) THEN
1.13 bertrand 636: SGLB = MAX( SGLB, TAU )
1.1 bertrand 637: ELSE
1.13 bertrand 638: SGUB = MIN( SGUB, TAU )
1.1 bertrand 639: END IF
640: *
641: * Calculate the new step
642: *
643: NITER = NITER + 1
644: IF( .NOT.SWTCH3 ) THEN
645: DTIPSQ = WORK( IP1 )*DELTA( IP1 )
646: DTISQ = WORK( I )*DELTA( I )
647: IF( ORGATI ) THEN
648: C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
649: ELSE
650: C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
651: END IF
652: A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
653: B = DTIPSQ*DTISQ*W
654: IF( C.EQ.ZERO ) THEN
655: IF( A.EQ.ZERO ) THEN
656: IF( ORGATI ) THEN
657: A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*( DPSI+DPHI )
658: ELSE
659: A = Z( IP1 )*Z( IP1 ) + DTISQ*DTISQ*( DPSI+DPHI )
660: END IF
661: END IF
662: ETA = B / A
663: ELSE IF( A.LE.ZERO ) THEN
664: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
665: ELSE
666: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
667: END IF
668: ELSE
669: *
670: * Interpolation using THREE most relevant poles
671: *
672: DTIIM = WORK( IIM1 )*DELTA( IIM1 )
673: DTIIP = WORK( IIP1 )*DELTA( IIP1 )
674: TEMP = RHOINV + PSI + PHI
675: IF( ORGATI ) THEN
676: TEMP1 = Z( IIM1 ) / DTIIM
677: TEMP1 = TEMP1*TEMP1
678: C = ( TEMP - DTIIP*( DPSI+DPHI ) ) -
679: $ ( D( IIM1 )-D( IIP1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
680: ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
681: IF( DPSI.LT.TEMP1 ) THEN
682: ZZ( 3 ) = DTIIP*DTIIP*DPHI
683: ELSE
684: ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
685: END IF
686: ELSE
687: TEMP1 = Z( IIP1 ) / DTIIP
688: TEMP1 = TEMP1*TEMP1
689: C = ( TEMP - DTIIM*( DPSI+DPHI ) ) -
690: $ ( D( IIP1 )-D( IIM1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
691: IF( DPHI.LT.TEMP1 ) THEN
692: ZZ( 1 ) = DTIIM*DTIIM*DPSI
693: ELSE
694: ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
695: END IF
696: ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
697: END IF
698: ZZ( 2 ) = Z( II )*Z( II )
699: DD( 1 ) = DTIIM
700: DD( 2 ) = DELTA( II )*WORK( II )
701: DD( 3 ) = DTIIP
702: CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
1.13 bertrand 703: *
704: IF( INFO.NE.0 ) THEN
705: *
706: * If INFO is not 0, i.e., DLAED6 failed, switch back
707: * to 2 pole interpolation.
708: *
709: SWTCH3 = .FALSE.
710: INFO = 0
711: DTIPSQ = WORK( IP1 )*DELTA( IP1 )
712: DTISQ = WORK( I )*DELTA( I )
713: IF( ORGATI ) THEN
714: C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
715: ELSE
716: C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
717: END IF
718: A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
719: B = DTIPSQ*DTISQ*W
720: IF( C.EQ.ZERO ) THEN
721: IF( A.EQ.ZERO ) THEN
722: IF( ORGATI ) THEN
723: A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*( DPSI+DPHI )
724: ELSE
725: A = Z( IP1 )*Z( IP1 ) + DTISQ*DTISQ*( DPSI+DPHI)
726: END IF
727: END IF
728: ETA = B / A
729: ELSE IF( A.LE.ZERO ) THEN
730: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
731: ELSE
732: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
733: END IF
734: END IF
1.1 bertrand 735: END IF
736: *
737: * Note, eta should be positive if w is negative, and
738: * eta should be negative otherwise. However,
739: * if for some reason caused by roundoff, eta*w > 0,
740: * we simply use one Newton step instead. This way
741: * will guarantee eta*w < 0.
742: *
743: IF( W*ETA.GE.ZERO )
744: $ ETA = -W / DW
1.13 bertrand 745: *
746: ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
747: TEMP = TAU + ETA
748: IF( TEMP.GT.SGUB .OR. TEMP.LT.SGLB ) THEN
1.1 bertrand 749: IF( W.LT.ZERO ) THEN
1.13 bertrand 750: ETA = ( SGUB-TAU ) / TWO
1.1 bertrand 751: ELSE
1.13 bertrand 752: ETA = ( SGLB-TAU ) / TWO
753: END IF
754: IF( GEOMAVG ) THEN
755: IF( W .LT. ZERO ) THEN
756: IF( TAU .GT. ZERO ) THEN
757: ETA = SQRT(SGUB*TAU)-TAU
758: END IF
759: ELSE
760: IF( SGLB .GT. ZERO ) THEN
761: ETA = SQRT(SGLB*TAU)-TAU
762: END IF
763: END IF
1.1 bertrand 764: END IF
765: END IF
766: *
767: PREW = W
768: *
1.13 bertrand 769: TAU = TAU + ETA
1.1 bertrand 770: SIGMA = SIGMA + ETA
1.13 bertrand 771: *
1.1 bertrand 772: DO 170 J = 1, N
773: WORK( J ) = WORK( J ) + ETA
774: DELTA( J ) = DELTA( J ) - ETA
775: 170 CONTINUE
776: *
777: * Evaluate PSI and the derivative DPSI
778: *
779: DPSI = ZERO
780: PSI = ZERO
781: ERRETM = ZERO
782: DO 180 J = 1, IIM1
783: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
784: PSI = PSI + Z( J )*TEMP
785: DPSI = DPSI + TEMP*TEMP
786: ERRETM = ERRETM + PSI
787: 180 CONTINUE
788: ERRETM = ABS( ERRETM )
789: *
790: * Evaluate PHI and the derivative DPHI
791: *
792: DPHI = ZERO
793: PHI = ZERO
794: DO 190 J = N, IIP1, -1
795: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
796: PHI = PHI + Z( J )*TEMP
797: DPHI = DPHI + TEMP*TEMP
798: ERRETM = ERRETM + PHI
799: 190 CONTINUE
800: *
1.13 bertrand 801: TAU2 = WORK( II )*DELTA( II )
802: TEMP = Z( II ) / TAU2
1.1 bertrand 803: DW = DPSI + DPHI + TEMP*TEMP
804: TEMP = Z( II )*TEMP
805: W = RHOINV + PHI + PSI + TEMP
1.13 bertrand 806: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
807: $ + THREE*ABS( TEMP )
808: * $ + ABS( TAU2 )*DW
1.1 bertrand 809: *
810: SWTCH = .FALSE.
811: IF( ORGATI ) THEN
812: IF( -W.GT.ABS( PREW ) / TEN )
813: $ SWTCH = .TRUE.
814: ELSE
815: IF( W.GT.ABS( PREW ) / TEN )
816: $ SWTCH = .TRUE.
817: END IF
818: *
819: * Main loop to update the values of the array DELTA and WORK
820: *
821: ITER = NITER + 1
822: *
823: DO 230 NITER = ITER, MAXIT
824: *
825: * Test for convergence
826: *
827: IF( ABS( W ).LE.EPS*ERRETM ) THEN
1.13 bertrand 828: * $ .OR. (SGUB-SGLB).LE.EIGHT*ABS(SGUB+SGLB) ) THEN
1.1 bertrand 829: GO TO 240
830: END IF
831: *
1.13 bertrand 832: IF( W.LE.ZERO ) THEN
833: SGLB = MAX( SGLB, TAU )
834: ELSE
835: SGUB = MIN( SGUB, TAU )
836: END IF
837: *
1.1 bertrand 838: * Calculate the new step
839: *
840: IF( .NOT.SWTCH3 ) THEN
841: DTIPSQ = WORK( IP1 )*DELTA( IP1 )
842: DTISQ = WORK( I )*DELTA( I )
843: IF( .NOT.SWTCH ) THEN
844: IF( ORGATI ) THEN
845: C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
846: ELSE
847: C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
848: END IF
849: ELSE
850: TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
851: IF( ORGATI ) THEN
852: DPSI = DPSI + TEMP*TEMP
853: ELSE
854: DPHI = DPHI + TEMP*TEMP
855: END IF
856: C = W - DTISQ*DPSI - DTIPSQ*DPHI
857: END IF
858: A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
859: B = DTIPSQ*DTISQ*W
860: IF( C.EQ.ZERO ) THEN
861: IF( A.EQ.ZERO ) THEN
862: IF( .NOT.SWTCH ) THEN
863: IF( ORGATI ) THEN
864: A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*
865: $ ( DPSI+DPHI )
866: ELSE
867: A = Z( IP1 )*Z( IP1 ) +
868: $ DTISQ*DTISQ*( DPSI+DPHI )
869: END IF
870: ELSE
871: A = DTISQ*DTISQ*DPSI + DTIPSQ*DTIPSQ*DPHI
872: END IF
873: END IF
874: ETA = B / A
875: ELSE IF( A.LE.ZERO ) THEN
876: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
877: ELSE
878: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
879: END IF
880: ELSE
881: *
882: * Interpolation using THREE most relevant poles
883: *
884: DTIIM = WORK( IIM1 )*DELTA( IIM1 )
885: DTIIP = WORK( IIP1 )*DELTA( IIP1 )
886: TEMP = RHOINV + PSI + PHI
887: IF( SWTCH ) THEN
888: C = TEMP - DTIIM*DPSI - DTIIP*DPHI
889: ZZ( 1 ) = DTIIM*DTIIM*DPSI
890: ZZ( 3 ) = DTIIP*DTIIP*DPHI
891: ELSE
892: IF( ORGATI ) THEN
893: TEMP1 = Z( IIM1 ) / DTIIM
894: TEMP1 = TEMP1*TEMP1
895: TEMP2 = ( D( IIM1 )-D( IIP1 ) )*
896: $ ( D( IIM1 )+D( IIP1 ) )*TEMP1
897: C = TEMP - DTIIP*( DPSI+DPHI ) - TEMP2
898: ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
899: IF( DPSI.LT.TEMP1 ) THEN
900: ZZ( 3 ) = DTIIP*DTIIP*DPHI
901: ELSE
902: ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
903: END IF
904: ELSE
905: TEMP1 = Z( IIP1 ) / DTIIP
906: TEMP1 = TEMP1*TEMP1
907: TEMP2 = ( D( IIP1 )-D( IIM1 ) )*
908: $ ( D( IIM1 )+D( IIP1 ) )*TEMP1
909: C = TEMP - DTIIM*( DPSI+DPHI ) - TEMP2
910: IF( DPHI.LT.TEMP1 ) THEN
911: ZZ( 1 ) = DTIIM*DTIIM*DPSI
912: ELSE
913: ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
914: END IF
915: ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
916: END IF
917: END IF
918: DD( 1 ) = DTIIM
919: DD( 2 ) = DELTA( II )*WORK( II )
920: DD( 3 ) = DTIIP
921: CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
1.13 bertrand 922: *
923: IF( INFO.NE.0 ) THEN
924: *
925: * If INFO is not 0, i.e., DLAED6 failed, switch
926: * back to two pole interpolation
927: *
928: SWTCH3 = .FALSE.
929: INFO = 0
930: DTIPSQ = WORK( IP1 )*DELTA( IP1 )
931: DTISQ = WORK( I )*DELTA( I )
932: IF( .NOT.SWTCH ) THEN
933: IF( ORGATI ) THEN
934: C = W - DTIPSQ*DW + DELSQ*( Z( I )/DTISQ )**2
935: ELSE
936: C = W - DTISQ*DW - DELSQ*( Z( IP1 )/DTIPSQ )**2
937: END IF
938: ELSE
939: TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
940: IF( ORGATI ) THEN
941: DPSI = DPSI + TEMP*TEMP
942: ELSE
943: DPHI = DPHI + TEMP*TEMP
944: END IF
945: C = W - DTISQ*DPSI - DTIPSQ*DPHI
946: END IF
947: A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
948: B = DTIPSQ*DTISQ*W
949: IF( C.EQ.ZERO ) THEN
950: IF( A.EQ.ZERO ) THEN
951: IF( .NOT.SWTCH ) THEN
952: IF( ORGATI ) THEN
953: A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*
954: $ ( DPSI+DPHI )
955: ELSE
956: A = Z( IP1 )*Z( IP1 ) +
957: $ DTISQ*DTISQ*( DPSI+DPHI )
958: END IF
959: ELSE
960: A = DTISQ*DTISQ*DPSI + DTIPSQ*DTIPSQ*DPHI
961: END IF
962: END IF
963: ETA = B / A
964: ELSE IF( A.LE.ZERO ) THEN
965: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
966: ELSE
967: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
968: END IF
969: END IF
1.1 bertrand 970: END IF
971: *
972: * Note, eta should be positive if w is negative, and
973: * eta should be negative otherwise. However,
974: * if for some reason caused by roundoff, eta*w > 0,
975: * we simply use one Newton step instead. This way
976: * will guarantee eta*w < 0.
977: *
978: IF( W*ETA.GE.ZERO )
979: $ ETA = -W / DW
1.13 bertrand 980: *
981: ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
982: TEMP=TAU+ETA
983: IF( TEMP.GT.SGUB .OR. TEMP.LT.SGLB ) THEN
1.1 bertrand 984: IF( W.LT.ZERO ) THEN
1.13 bertrand 985: ETA = ( SGUB-TAU ) / TWO
1.1 bertrand 986: ELSE
1.13 bertrand 987: ETA = ( SGLB-TAU ) / TWO
988: END IF
989: IF( GEOMAVG ) THEN
990: IF( W .LT. ZERO ) THEN
991: IF( TAU .GT. ZERO ) THEN
992: ETA = SQRT(SGUB*TAU)-TAU
993: END IF
994: ELSE
995: IF( SGLB .GT. ZERO ) THEN
996: ETA = SQRT(SGLB*TAU)-TAU
997: END IF
998: END IF
1.1 bertrand 999: END IF
1000: END IF
1001: *
1.13 bertrand 1002: PREW = W
1003: *
1.1 bertrand 1004: TAU = TAU + ETA
1.13 bertrand 1005: SIGMA = SIGMA + ETA
1.1 bertrand 1006: *
1007: DO 200 J = 1, N
1008: WORK( J ) = WORK( J ) + ETA
1009: DELTA( J ) = DELTA( J ) - ETA
1010: 200 CONTINUE
1011: *
1012: * Evaluate PSI and the derivative DPSI
1013: *
1014: DPSI = ZERO
1015: PSI = ZERO
1016: ERRETM = ZERO
1017: DO 210 J = 1, IIM1
1018: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
1019: PSI = PSI + Z( J )*TEMP
1020: DPSI = DPSI + TEMP*TEMP
1021: ERRETM = ERRETM + PSI
1022: 210 CONTINUE
1023: ERRETM = ABS( ERRETM )
1024: *
1025: * Evaluate PHI and the derivative DPHI
1026: *
1027: DPHI = ZERO
1028: PHI = ZERO
1029: DO 220 J = N, IIP1, -1
1030: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
1031: PHI = PHI + Z( J )*TEMP
1032: DPHI = DPHI + TEMP*TEMP
1033: ERRETM = ERRETM + PHI
1034: 220 CONTINUE
1035: *
1.13 bertrand 1036: TAU2 = WORK( II )*DELTA( II )
1037: TEMP = Z( II ) / TAU2
1.1 bertrand 1038: DW = DPSI + DPHI + TEMP*TEMP
1039: TEMP = Z( II )*TEMP
1040: W = RHOINV + PHI + PSI + TEMP
1.13 bertrand 1041: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
1042: $ + THREE*ABS( TEMP )
1043: * $ + ABS( TAU2 )*DW
1044: *
1.1 bertrand 1045: IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )
1046: $ SWTCH = .NOT.SWTCH
1047: *
1048: 230 CONTINUE
1049: *
1050: * Return with INFO = 1, NITER = MAXIT and not converged
1051: *
1052: INFO = 1
1053: *
1054: END IF
1055: *
1056: 240 CONTINUE
1057: RETURN
1058: *
1059: * End of DLASD4
1060: *
1061: 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