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