1: *> \brief \b DLASD4
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 November 2011
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.0) --
157: * -- LAPACK is a software package provided by Univ. of Tennessee, --
158: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159: * November 2011
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 = 64 )
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
181: INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
182: DOUBLE PRECISION A, B, C, DELSQ, DELSQ2, DPHI, DPSI, DTIIM,
183: $ DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS,
184: $ ERRETM, ETA, PHI, PREW, PSI, RHOINV, SG2LB,
185: $ SG2UB, TAU, 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 TAU 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: TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
275: ELSE
276: TAU = ( 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+TAU <= 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 TAU is to approximate
289: * SIGMA_n^2 - D( N )*D( N )
290: *
291: IF( A.LT.ZERO ) THEN
292: TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
293: ELSE
294: TAU = ( 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+TAU < SIGMA(N)^2 < D(N)^2+RHO/2
299: *
300: END IF
301: *
302: * The following ETA is to approximate SIGMA_n - D( N )
303: *
304: ETA = TAU / ( D( N )+SQRT( D( N )*D( N )+TAU ) )
305: *
306: SIGMA = D( N ) + ETA
307: DO 30 J = 1, N
308: DELTA( J ) = ( D( J )-D( I ) ) - ETA
309: WORK( J ) = D( J ) + D( I ) + ETA
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( TAU )*( 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: TAU = TAU + ETA
372: ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
373: DO 50 J = 1, N
374: DELTA( J ) = DELTA( J ) - ETA
375: WORK( J ) = WORK( J ) + ETA
376: 50 CONTINUE
377: *
378: SIGMA = SIGMA + ETA
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: TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )
396: PHI = Z( N )*TEMP
397: DPHI = TEMP*TEMP
398: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
399: $ ABS( TAU )*( DPSI+DPHI )
400: *
401: W = RHOINV + PHI + PSI
402: *
403: * Main loop to update the values of the array DELTA
404: *
405: ITER = NITER + 1
406: *
407: DO 90 NITER = ITER, MAXIT
408: *
409: * Test for convergence
410: *
411: IF( ABS( W ).LE.EPS*ERRETM ) THEN
412: GO TO 240
413: END IF
414: *
415: * Calculate the new step
416: *
417: DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
418: DTNSQ = WORK( N )*DELTA( N )
419: C = W - DTNSQ1*DPSI - DTNSQ*DPHI
420: A = ( DTNSQ+DTNSQ1 )*W - DTNSQ1*DTNSQ*( DPSI+DPHI )
421: B = DTNSQ1*DTNSQ*W
422: IF( A.GE.ZERO ) THEN
423: ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
424: ELSE
425: ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
426: END IF
427: *
428: * Note, eta should be positive if w is negative, and
429: * eta should be negative otherwise. However,
430: * if for some reason caused by roundoff, eta*w > 0,
431: * we simply use one Newton step instead. This way
432: * will guarantee eta*w < 0.
433: *
434: IF( W*ETA.GT.ZERO )
435: $ ETA = -W / ( DPSI+DPHI )
436: TEMP = ETA - DTNSQ
437: IF( TEMP.LE.ZERO )
438: $ ETA = ETA / TWO
439: *
440: TAU = TAU + ETA
441: ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
442: DO 70 J = 1, N
443: DELTA( J ) = DELTA( J ) - ETA
444: WORK( J ) = WORK( J ) + ETA
445: 70 CONTINUE
446: *
447: SIGMA = SIGMA + ETA
448: *
449: * Evaluate PSI and the derivative DPSI
450: *
451: DPSI = ZERO
452: PSI = ZERO
453: ERRETM = ZERO
454: DO 80 J = 1, II
455: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
456: PSI = PSI + Z( J )*TEMP
457: DPSI = DPSI + TEMP*TEMP
458: ERRETM = ERRETM + PSI
459: 80 CONTINUE
460: ERRETM = ABS( ERRETM )
461: *
462: * Evaluate PHI and the derivative DPHI
463: *
464: TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )
465: PHI = Z( N )*TEMP
466: DPHI = TEMP*TEMP
467: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
468: $ ABS( TAU )*( DPSI+DPHI )
469: *
470: W = RHOINV + PHI + PSI
471: 90 CONTINUE
472: *
473: * Return with INFO = 1, NITER = MAXIT and not converged
474: *
475: INFO = 1
476: GO TO 240
477: *
478: * End for the case I = N
479: *
480: ELSE
481: *
482: * The case for I < N
483: *
484: NITER = 1
485: IP1 = I + 1
486: *
487: * Calculate initial guess
488: *
489: DELSQ = ( D( IP1 )-D( I ) )*( D( IP1 )+D( I ) )
490: DELSQ2 = DELSQ / TWO
491: TEMP = DELSQ2 / ( D( I )+SQRT( D( I )*D( I )+DELSQ2 ) )
492: DO 100 J = 1, N
493: WORK( J ) = D( J ) + D( I ) + TEMP
494: DELTA( J ) = ( D( J )-D( I ) ) - TEMP
495: 100 CONTINUE
496: *
497: PSI = ZERO
498: DO 110 J = 1, I - 1
499: PSI = PSI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
500: 110 CONTINUE
501: *
502: PHI = ZERO
503: DO 120 J = N, I + 2, -1
504: PHI = PHI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
505: 120 CONTINUE
506: C = RHOINV + PSI + PHI
507: W = C + Z( I )*Z( I ) / ( WORK( I )*DELTA( I ) ) +
508: $ Z( IP1 )*Z( IP1 ) / ( WORK( IP1 )*DELTA( IP1 ) )
509: *
510: IF( W.GT.ZERO ) THEN
511: *
512: * d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2
513: *
514: * We choose d(i) as origin.
515: *
516: ORGATI = .TRUE.
517: SG2LB = ZERO
518: SG2UB = DELSQ2
519: A = C*DELSQ + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
520: B = Z( I )*Z( I )*DELSQ
521: IF( A.GT.ZERO ) THEN
522: TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
523: ELSE
524: TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
525: END IF
526: *
527: * TAU now is an estimation of SIGMA^2 - D( I )^2. The
528: * following, however, is the corresponding estimation of
529: * SIGMA - D( I ).
530: *
531: ETA = TAU / ( D( I )+SQRT( D( I )*D( I )+TAU ) )
532: ELSE
533: *
534: * (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2
535: *
536: * We choose d(i+1) as origin.
537: *
538: ORGATI = .FALSE.
539: SG2LB = -DELSQ2
540: SG2UB = ZERO
541: A = C*DELSQ - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
542: B = Z( IP1 )*Z( IP1 )*DELSQ
543: IF( A.LT.ZERO ) THEN
544: TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
545: ELSE
546: TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
547: END IF
548: *
549: * TAU now is an estimation of SIGMA^2 - D( IP1 )^2. The
550: * following, however, is the corresponding estimation of
551: * SIGMA - D( IP1 ).
552: *
553: ETA = TAU / ( D( IP1 )+SQRT( ABS( D( IP1 )*D( IP1 )+
554: $ TAU ) ) )
555: END IF
556: *
557: IF( ORGATI ) THEN
558: II = I
559: SIGMA = D( I ) + ETA
560: DO 130 J = 1, N
561: WORK( J ) = D( J ) + D( I ) + ETA
562: DELTA( J ) = ( D( J )-D( I ) ) - ETA
563: 130 CONTINUE
564: ELSE
565: II = I + 1
566: SIGMA = D( IP1 ) + ETA
567: DO 140 J = 1, N
568: WORK( J ) = D( J ) + D( IP1 ) + ETA
569: DELTA( J ) = ( D( J )-D( IP1 ) ) - ETA
570: 140 CONTINUE
571: END IF
572: IIM1 = II - 1
573: IIP1 = II + 1
574: *
575: * Evaluate PSI and the derivative DPSI
576: *
577: DPSI = ZERO
578: PSI = ZERO
579: ERRETM = ZERO
580: DO 150 J = 1, IIM1
581: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
582: PSI = PSI + Z( J )*TEMP
583: DPSI = DPSI + TEMP*TEMP
584: ERRETM = ERRETM + PSI
585: 150 CONTINUE
586: ERRETM = ABS( ERRETM )
587: *
588: * Evaluate PHI and the derivative DPHI
589: *
590: DPHI = ZERO
591: PHI = ZERO
592: DO 160 J = N, IIP1, -1
593: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
594: PHI = PHI + Z( J )*TEMP
595: DPHI = DPHI + TEMP*TEMP
596: ERRETM = ERRETM + PHI
597: 160 CONTINUE
598: *
599: W = RHOINV + PHI + PSI
600: *
601: * W is the value of the secular function with
602: * its ii-th element removed.
603: *
604: SWTCH3 = .FALSE.
605: IF( ORGATI ) THEN
606: IF( W.LT.ZERO )
607: $ SWTCH3 = .TRUE.
608: ELSE
609: IF( W.GT.ZERO )
610: $ SWTCH3 = .TRUE.
611: END IF
612: IF( II.EQ.1 .OR. II.EQ.N )
613: $ SWTCH3 = .FALSE.
614: *
615: TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
616: DW = DPSI + DPHI + TEMP*TEMP
617: TEMP = Z( II )*TEMP
618: W = W + TEMP
619: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
620: $ THREE*ABS( TEMP ) + ABS( TAU )*DW
621: *
622: * Test for convergence
623: *
624: IF( ABS( W ).LE.EPS*ERRETM ) THEN
625: GO TO 240
626: END IF
627: *
628: IF( W.LE.ZERO ) THEN
629: SG2LB = MAX( SG2LB, TAU )
630: ELSE
631: SG2UB = MIN( SG2UB, TAU )
632: END IF
633: *
634: * Calculate the new step
635: *
636: NITER = NITER + 1
637: IF( .NOT.SWTCH3 ) THEN
638: DTIPSQ = WORK( IP1 )*DELTA( IP1 )
639: DTISQ = WORK( I )*DELTA( I )
640: IF( ORGATI ) THEN
641: C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
642: ELSE
643: C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
644: END IF
645: A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
646: B = DTIPSQ*DTISQ*W
647: IF( C.EQ.ZERO ) THEN
648: IF( A.EQ.ZERO ) THEN
649: IF( ORGATI ) THEN
650: A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*( DPSI+DPHI )
651: ELSE
652: A = Z( IP1 )*Z( IP1 ) + DTISQ*DTISQ*( DPSI+DPHI )
653: END IF
654: END IF
655: ETA = B / A
656: ELSE IF( A.LE.ZERO ) THEN
657: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
658: ELSE
659: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
660: END IF
661: ELSE
662: *
663: * Interpolation using THREE most relevant poles
664: *
665: DTIIM = WORK( IIM1 )*DELTA( IIM1 )
666: DTIIP = WORK( IIP1 )*DELTA( IIP1 )
667: TEMP = RHOINV + PSI + PHI
668: IF( ORGATI ) THEN
669: TEMP1 = Z( IIM1 ) / DTIIM
670: TEMP1 = TEMP1*TEMP1
671: C = ( TEMP - DTIIP*( DPSI+DPHI ) ) -
672: $ ( D( IIM1 )-D( IIP1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
673: ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
674: IF( DPSI.LT.TEMP1 ) THEN
675: ZZ( 3 ) = DTIIP*DTIIP*DPHI
676: ELSE
677: ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
678: END IF
679: ELSE
680: TEMP1 = Z( IIP1 ) / DTIIP
681: TEMP1 = TEMP1*TEMP1
682: C = ( TEMP - DTIIM*( DPSI+DPHI ) ) -
683: $ ( D( IIP1 )-D( IIM1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
684: IF( DPHI.LT.TEMP1 ) THEN
685: ZZ( 1 ) = DTIIM*DTIIM*DPSI
686: ELSE
687: ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
688: END IF
689: ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
690: END IF
691: ZZ( 2 ) = Z( II )*Z( II )
692: DD( 1 ) = DTIIM
693: DD( 2 ) = DELTA( II )*WORK( II )
694: DD( 3 ) = DTIIP
695: CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
696: IF( INFO.NE.0 )
697: $ GO TO 240
698: END IF
699: *
700: * Note, eta should be positive if w is negative, and
701: * eta should be negative otherwise. However,
702: * if for some reason caused by roundoff, eta*w > 0,
703: * we simply use one Newton step instead. This way
704: * will guarantee eta*w < 0.
705: *
706: IF( W*ETA.GE.ZERO )
707: $ ETA = -W / DW
708: IF( ORGATI ) THEN
709: TEMP1 = WORK( I )*DELTA( I )
710: TEMP = ETA - TEMP1
711: ELSE
712: TEMP1 = WORK( IP1 )*DELTA( IP1 )
713: TEMP = ETA - TEMP1
714: END IF
715: IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN
716: IF( W.LT.ZERO ) THEN
717: ETA = ( SG2UB-TAU ) / TWO
718: ELSE
719: ETA = ( SG2LB-TAU ) / TWO
720: END IF
721: END IF
722: *
723: TAU = TAU + ETA
724: ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
725: *
726: PREW = W
727: *
728: SIGMA = SIGMA + ETA
729: DO 170 J = 1, N
730: WORK( J ) = WORK( J ) + ETA
731: DELTA( J ) = DELTA( J ) - ETA
732: 170 CONTINUE
733: *
734: * Evaluate PSI and the derivative DPSI
735: *
736: DPSI = ZERO
737: PSI = ZERO
738: ERRETM = ZERO
739: DO 180 J = 1, IIM1
740: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
741: PSI = PSI + Z( J )*TEMP
742: DPSI = DPSI + TEMP*TEMP
743: ERRETM = ERRETM + PSI
744: 180 CONTINUE
745: ERRETM = ABS( ERRETM )
746: *
747: * Evaluate PHI and the derivative DPHI
748: *
749: DPHI = ZERO
750: PHI = ZERO
751: DO 190 J = N, IIP1, -1
752: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
753: PHI = PHI + Z( J )*TEMP
754: DPHI = DPHI + TEMP*TEMP
755: ERRETM = ERRETM + PHI
756: 190 CONTINUE
757: *
758: TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
759: DW = DPSI + DPHI + TEMP*TEMP
760: TEMP = Z( II )*TEMP
761: W = RHOINV + PHI + PSI + TEMP
762: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
763: $ THREE*ABS( TEMP ) + ABS( TAU )*DW
764: *
765: IF( W.LE.ZERO ) THEN
766: SG2LB = MAX( SG2LB, TAU )
767: ELSE
768: SG2UB = MIN( SG2UB, TAU )
769: END IF
770: *
771: SWTCH = .FALSE.
772: IF( ORGATI ) THEN
773: IF( -W.GT.ABS( PREW ) / TEN )
774: $ SWTCH = .TRUE.
775: ELSE
776: IF( W.GT.ABS( PREW ) / TEN )
777: $ SWTCH = .TRUE.
778: END IF
779: *
780: * Main loop to update the values of the array DELTA and WORK
781: *
782: ITER = NITER + 1
783: *
784: DO 230 NITER = ITER, MAXIT
785: *
786: * Test for convergence
787: *
788: IF( ABS( W ).LE.EPS*ERRETM ) THEN
789: GO TO 240
790: END IF
791: *
792: * Calculate the new step
793: *
794: IF( .NOT.SWTCH3 ) THEN
795: DTIPSQ = WORK( IP1 )*DELTA( IP1 )
796: DTISQ = WORK( I )*DELTA( I )
797: IF( .NOT.SWTCH ) THEN
798: IF( ORGATI ) THEN
799: C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
800: ELSE
801: C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
802: END IF
803: ELSE
804: TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
805: IF( ORGATI ) THEN
806: DPSI = DPSI + TEMP*TEMP
807: ELSE
808: DPHI = DPHI + TEMP*TEMP
809: END IF
810: C = W - DTISQ*DPSI - DTIPSQ*DPHI
811: END IF
812: A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
813: B = DTIPSQ*DTISQ*W
814: IF( C.EQ.ZERO ) THEN
815: IF( A.EQ.ZERO ) THEN
816: IF( .NOT.SWTCH ) THEN
817: IF( ORGATI ) THEN
818: A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*
819: $ ( DPSI+DPHI )
820: ELSE
821: A = Z( IP1 )*Z( IP1 ) +
822: $ DTISQ*DTISQ*( DPSI+DPHI )
823: END IF
824: ELSE
825: A = DTISQ*DTISQ*DPSI + DTIPSQ*DTIPSQ*DPHI
826: END IF
827: END IF
828: ETA = B / A
829: ELSE IF( A.LE.ZERO ) THEN
830: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
831: ELSE
832: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
833: END IF
834: ELSE
835: *
836: * Interpolation using THREE most relevant poles
837: *
838: DTIIM = WORK( IIM1 )*DELTA( IIM1 )
839: DTIIP = WORK( IIP1 )*DELTA( IIP1 )
840: TEMP = RHOINV + PSI + PHI
841: IF( SWTCH ) THEN
842: C = TEMP - DTIIM*DPSI - DTIIP*DPHI
843: ZZ( 1 ) = DTIIM*DTIIM*DPSI
844: ZZ( 3 ) = DTIIP*DTIIP*DPHI
845: ELSE
846: IF( ORGATI ) THEN
847: TEMP1 = Z( IIM1 ) / DTIIM
848: TEMP1 = TEMP1*TEMP1
849: TEMP2 = ( D( IIM1 )-D( IIP1 ) )*
850: $ ( D( IIM1 )+D( IIP1 ) )*TEMP1
851: C = TEMP - DTIIP*( DPSI+DPHI ) - TEMP2
852: ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
853: IF( DPSI.LT.TEMP1 ) THEN
854: ZZ( 3 ) = DTIIP*DTIIP*DPHI
855: ELSE
856: ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
857: END IF
858: ELSE
859: TEMP1 = Z( IIP1 ) / DTIIP
860: TEMP1 = TEMP1*TEMP1
861: TEMP2 = ( D( IIP1 )-D( IIM1 ) )*
862: $ ( D( IIM1 )+D( IIP1 ) )*TEMP1
863: C = TEMP - DTIIM*( DPSI+DPHI ) - TEMP2
864: IF( DPHI.LT.TEMP1 ) THEN
865: ZZ( 1 ) = DTIIM*DTIIM*DPSI
866: ELSE
867: ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
868: END IF
869: ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
870: END IF
871: END IF
872: DD( 1 ) = DTIIM
873: DD( 2 ) = DELTA( II )*WORK( II )
874: DD( 3 ) = DTIIP
875: CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
876: IF( INFO.NE.0 )
877: $ GO TO 240
878: END IF
879: *
880: * Note, eta should be positive if w is negative, and
881: * eta should be negative otherwise. However,
882: * if for some reason caused by roundoff, eta*w > 0,
883: * we simply use one Newton step instead. This way
884: * will guarantee eta*w < 0.
885: *
886: IF( W*ETA.GE.ZERO )
887: $ ETA = -W / DW
888: IF( ORGATI ) THEN
889: TEMP1 = WORK( I )*DELTA( I )
890: TEMP = ETA - TEMP1
891: ELSE
892: TEMP1 = WORK( IP1 )*DELTA( IP1 )
893: TEMP = ETA - TEMP1
894: END IF
895: IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN
896: IF( W.LT.ZERO ) THEN
897: ETA = ( SG2UB-TAU ) / TWO
898: ELSE
899: ETA = ( SG2LB-TAU ) / TWO
900: END IF
901: END IF
902: *
903: TAU = TAU + ETA
904: ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
905: *
906: SIGMA = SIGMA + ETA
907: DO 200 J = 1, N
908: WORK( J ) = WORK( J ) + ETA
909: DELTA( J ) = DELTA( J ) - ETA
910: 200 CONTINUE
911: *
912: PREW = W
913: *
914: * Evaluate PSI and the derivative DPSI
915: *
916: DPSI = ZERO
917: PSI = ZERO
918: ERRETM = ZERO
919: DO 210 J = 1, IIM1
920: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
921: PSI = PSI + Z( J )*TEMP
922: DPSI = DPSI + TEMP*TEMP
923: ERRETM = ERRETM + PSI
924: 210 CONTINUE
925: ERRETM = ABS( ERRETM )
926: *
927: * Evaluate PHI and the derivative DPHI
928: *
929: DPHI = ZERO
930: PHI = ZERO
931: DO 220 J = N, IIP1, -1
932: TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
933: PHI = PHI + Z( J )*TEMP
934: DPHI = DPHI + TEMP*TEMP
935: ERRETM = ERRETM + PHI
936: 220 CONTINUE
937: *
938: TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
939: DW = DPSI + DPHI + TEMP*TEMP
940: TEMP = Z( II )*TEMP
941: W = RHOINV + PHI + PSI + TEMP
942: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
943: $ THREE*ABS( TEMP ) + ABS( TAU )*DW
944: IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )
945: $ SWTCH = .NOT.SWTCH
946: *
947: IF( W.LE.ZERO ) THEN
948: SG2LB = MAX( SG2LB, TAU )
949: ELSE
950: SG2UB = MIN( SG2UB, TAU )
951: END IF
952: *
953: 230 CONTINUE
954: *
955: * Return with INFO = 1, NITER = MAXIT and not converged
956: *
957: INFO = 1
958: *
959: END IF
960: *
961: 240 CONTINUE
962: RETURN
963: *
964: * End of DLASD4
965: *
966: 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