1: *> \brief \b DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAED6 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed6.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed6.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed6.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
22: *
23: * .. Scalar Arguments ..
24: * LOGICAL ORGATI
25: * INTEGER INFO, KNITER
26: * DOUBLE PRECISION FINIT, RHO, TAU
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION D( 3 ), Z( 3 )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLAED6 computes the positive or negative root (closest to the origin)
39: *> of
40: *> z(1) z(2) z(3)
41: *> f(x) = rho + --------- + ---------- + ---------
42: *> d(1)-x d(2)-x d(3)-x
43: *>
44: *> It is assumed that
45: *>
46: *> if ORGATI = .true. the root is between d(2) and d(3);
47: *> otherwise it is between d(1) and d(2)
48: *>
49: *> This routine will be called by DLAED4 when necessary. In most cases,
50: *> the root sought is the smallest in magnitude, though it might not be
51: *> in some extremely rare situations.
52: *> \endverbatim
53: *
54: * Arguments:
55: * ==========
56: *
57: *> \param[in] KNITER
58: *> \verbatim
59: *> KNITER is INTEGER
60: *> Refer to DLAED4 for its significance.
61: *> \endverbatim
62: *>
63: *> \param[in] ORGATI
64: *> \verbatim
65: *> ORGATI is LOGICAL
66: *> If ORGATI is true, the needed root is between d(2) and
67: *> d(3); otherwise it is between d(1) and d(2). See
68: *> DLAED4 for further details.
69: *> \endverbatim
70: *>
71: *> \param[in] RHO
72: *> \verbatim
73: *> RHO is DOUBLE PRECISION
74: *> Refer to the equation f(x) above.
75: *> \endverbatim
76: *>
77: *> \param[in] D
78: *> \verbatim
79: *> D is DOUBLE PRECISION array, dimension (3)
80: *> D satisfies d(1) < d(2) < d(3).
81: *> \endverbatim
82: *>
83: *> \param[in] Z
84: *> \verbatim
85: *> Z is DOUBLE PRECISION array, dimension (3)
86: *> Each of the elements in z must be positive.
87: *> \endverbatim
88: *>
89: *> \param[in] FINIT
90: *> \verbatim
91: *> FINIT is DOUBLE PRECISION
92: *> The value of f at 0. It is more accurate than the one
93: *> evaluated inside this routine (if someone wants to do
94: *> so).
95: *> \endverbatim
96: *>
97: *> \param[out] TAU
98: *> \verbatim
99: *> TAU is DOUBLE PRECISION
100: *> The root of the equation f(x).
101: *> \endverbatim
102: *>
103: *> \param[out] INFO
104: *> \verbatim
105: *> INFO is INTEGER
106: *> = 0: successful exit
107: *> > 0: if INFO = 1, failure to converge
108: *> \endverbatim
109: *
110: * Authors:
111: * ========
112: *
113: *> \author Univ. of Tennessee
114: *> \author Univ. of California Berkeley
115: *> \author Univ. of Colorado Denver
116: *> \author NAG Ltd.
117: *
118: *> \ingroup auxOTHERcomputational
119: *
120: *> \par Further Details:
121: * =====================
122: *>
123: *> \verbatim
124: *>
125: *> 10/02/03: This version has a few statements commented out for thread
126: *> safety (machine parameters are computed on each entry). SJH.
127: *>
128: *> 05/10/06: Modified from a new version of Ren-Cang Li, use
129: *> Gragg-Thornton-Warner cubic convergent scheme for better stability.
130: *> \endverbatim
131: *
132: *> \par Contributors:
133: * ==================
134: *>
135: *> Ren-Cang Li, Computer Science Division, University of California
136: *> at Berkeley, USA
137: *>
138: * =====================================================================
139: SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
140: *
141: * -- LAPACK computational routine --
142: * -- LAPACK is a software package provided by Univ. of Tennessee, --
143: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144: *
145: * .. Scalar Arguments ..
146: LOGICAL ORGATI
147: INTEGER INFO, KNITER
148: DOUBLE PRECISION FINIT, RHO, TAU
149: * ..
150: * .. Array Arguments ..
151: DOUBLE PRECISION D( 3 ), Z( 3 )
152: * ..
153: *
154: * =====================================================================
155: *
156: * .. Parameters ..
157: INTEGER MAXIT
158: PARAMETER ( MAXIT = 40 )
159: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT
160: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
161: $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 )
162: * ..
163: * .. External Functions ..
164: DOUBLE PRECISION DLAMCH
165: EXTERNAL DLAMCH
166: * ..
167: * .. Local Arrays ..
168: DOUBLE PRECISION DSCALE( 3 ), ZSCALE( 3 )
169: * ..
170: * .. Local Scalars ..
171: LOGICAL SCALE
172: INTEGER I, ITER, NITER
173: DOUBLE PRECISION A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
174: $ FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
175: $ SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4,
176: $ LBD, UBD
177: * ..
178: * .. Intrinsic Functions ..
179: INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT
180: * ..
181: * .. Executable Statements ..
182: *
183: INFO = 0
184: *
185: IF( ORGATI ) THEN
186: LBD = D(2)
187: UBD = D(3)
188: ELSE
189: LBD = D(1)
190: UBD = D(2)
191: END IF
192: IF( FINIT .LT. ZERO )THEN
193: LBD = ZERO
194: ELSE
195: UBD = ZERO
196: END IF
197: *
198: NITER = 1
199: TAU = ZERO
200: IF( KNITER.EQ.2 ) THEN
201: IF( ORGATI ) THEN
202: TEMP = ( D( 3 )-D( 2 ) ) / TWO
203: C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
204: A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
205: B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
206: ELSE
207: TEMP = ( D( 1 )-D( 2 ) ) / TWO
208: C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
209: A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
210: B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
211: END IF
212: TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
213: A = A / TEMP
214: B = B / TEMP
215: C = C / TEMP
216: IF( C.EQ.ZERO ) THEN
217: TAU = B / A
218: ELSE IF( A.LE.ZERO ) THEN
219: TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
220: ELSE
221: TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
222: END IF
223: IF( TAU .LT. LBD .OR. TAU .GT. UBD )
224: $ TAU = ( LBD+UBD )/TWO
225: IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN
226: TAU = ZERO
227: ELSE
228: TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) +
229: $ TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) +
230: $ TAU*Z(3)/( D(3)*( D( 3 )-TAU ) )
231: IF( TEMP .LE. ZERO )THEN
232: LBD = TAU
233: ELSE
234: UBD = TAU
235: END IF
236: IF( ABS( FINIT ).LE.ABS( TEMP ) )
237: $ TAU = ZERO
238: END IF
239: END IF
240: *
241: * get machine parameters for possible scaling to avoid overflow
242: *
243: * modified by Sven: parameters SMALL1, SMINV1, SMALL2,
244: * SMINV2, EPS are not SAVEd anymore between one call to the
245: * others but recomputed at each call
246: *
247: EPS = DLAMCH( 'Epsilon' )
248: BASE = DLAMCH( 'Base' )
249: SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) /
250: $ THREE ) )
251: SMINV1 = ONE / SMALL1
252: SMALL2 = SMALL1*SMALL1
253: SMINV2 = SMINV1*SMINV1
254: *
255: * Determine if scaling of inputs necessary to avoid overflow
256: * when computing 1/TEMP**3
257: *
258: IF( ORGATI ) THEN
259: TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
260: ELSE
261: TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
262: END IF
263: SCALE = .FALSE.
264: IF( TEMP.LE.SMALL1 ) THEN
265: SCALE = .TRUE.
266: IF( TEMP.LE.SMALL2 ) THEN
267: *
268: * Scale up by power of radix nearest 1/SAFMIN**(2/3)
269: *
270: SCLFAC = SMINV2
271: SCLINV = SMALL2
272: ELSE
273: *
274: * Scale up by power of radix nearest 1/SAFMIN**(1/3)
275: *
276: SCLFAC = SMINV1
277: SCLINV = SMALL1
278: END IF
279: *
280: * Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
281: *
282: DO 10 I = 1, 3
283: DSCALE( I ) = D( I )*SCLFAC
284: ZSCALE( I ) = Z( I )*SCLFAC
285: 10 CONTINUE
286: TAU = TAU*SCLFAC
287: LBD = LBD*SCLFAC
288: UBD = UBD*SCLFAC
289: ELSE
290: *
291: * Copy D and Z to DSCALE and ZSCALE
292: *
293: DO 20 I = 1, 3
294: DSCALE( I ) = D( I )
295: ZSCALE( I ) = Z( I )
296: 20 CONTINUE
297: END IF
298: *
299: FC = ZERO
300: DF = ZERO
301: DDF = ZERO
302: DO 30 I = 1, 3
303: TEMP = ONE / ( DSCALE( I )-TAU )
304: TEMP1 = ZSCALE( I )*TEMP
305: TEMP2 = TEMP1*TEMP
306: TEMP3 = TEMP2*TEMP
307: FC = FC + TEMP1 / DSCALE( I )
308: DF = DF + TEMP2
309: DDF = DDF + TEMP3
310: 30 CONTINUE
311: F = FINIT + TAU*FC
312: *
313: IF( ABS( F ).LE.ZERO )
314: $ GO TO 60
315: IF( F .LE. ZERO )THEN
316: LBD = TAU
317: ELSE
318: UBD = TAU
319: END IF
320: *
321: * Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
322: * scheme
323: *
324: * It is not hard to see that
325: *
326: * 1) Iterations will go up monotonically
327: * if FINIT < 0;
328: *
329: * 2) Iterations will go down monotonically
330: * if FINIT > 0.
331: *
332: ITER = NITER + 1
333: *
334: DO 50 NITER = ITER, MAXIT
335: *
336: IF( ORGATI ) THEN
337: TEMP1 = DSCALE( 2 ) - TAU
338: TEMP2 = DSCALE( 3 ) - TAU
339: ELSE
340: TEMP1 = DSCALE( 1 ) - TAU
341: TEMP2 = DSCALE( 2 ) - TAU
342: END IF
343: A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
344: B = TEMP1*TEMP2*F
345: C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
346: TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
347: A = A / TEMP
348: B = B / TEMP
349: C = C / TEMP
350: IF( C.EQ.ZERO ) THEN
351: ETA = B / A
352: ELSE IF( A.LE.ZERO ) THEN
353: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
354: ELSE
355: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
356: END IF
357: IF( F*ETA.GE.ZERO ) THEN
358: ETA = -F / DF
359: END IF
360: *
361: TAU = TAU + ETA
362: IF( TAU .LT. LBD .OR. TAU .GT. UBD )
363: $ TAU = ( LBD + UBD )/TWO
364: *
365: FC = ZERO
366: ERRETM = ZERO
367: DF = ZERO
368: DDF = ZERO
369: DO 40 I = 1, 3
370: IF ( ( DSCALE( I )-TAU ).NE.ZERO ) THEN
371: TEMP = ONE / ( DSCALE( I )-TAU )
372: TEMP1 = ZSCALE( I )*TEMP
373: TEMP2 = TEMP1*TEMP
374: TEMP3 = TEMP2*TEMP
375: TEMP4 = TEMP1 / DSCALE( I )
376: FC = FC + TEMP4
377: ERRETM = ERRETM + ABS( TEMP4 )
378: DF = DF + TEMP2
379: DDF = DDF + TEMP3
380: ELSE
381: GO TO 60
382: END IF
383: 40 CONTINUE
384: F = FINIT + TAU*FC
385: ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
386: $ ABS( TAU )*DF
387: IF( ( ABS( F ).LE.FOUR*EPS*ERRETM ) .OR.
388: $ ( (UBD-LBD).LE.FOUR*EPS*ABS(TAU) ) )
389: $ GO TO 60
390: IF( F .LE. ZERO )THEN
391: LBD = TAU
392: ELSE
393: UBD = TAU
394: END IF
395: 50 CONTINUE
396: INFO = 1
397: 60 CONTINUE
398: *
399: * Undo scaling
400: *
401: IF( SCALE )
402: $ TAU = TAU*SCLINV
403: RETURN
404: *
405: * End of DLAED6
406: *
407: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed6.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack