1: *> \brief \b DLAED6 used by sstedc. 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: *> \date November 2015
119: *
120: *> \ingroup auxOTHERcomputational
121: *
122: *> \par Further Details:
123: * =====================
124: *>
125: *> \verbatim
126: *>
127: *> 10/02/03: This version has a few statements commented out for thread
128: *> safety (machine parameters are computed on each entry). SJH.
129: *>
130: *> 05/10/06: Modified from a new version of Ren-Cang Li, use
131: *> Gragg-Thornton-Warner cubic convergent scheme for better stability.
132: *> \endverbatim
133: *
134: *> \par Contributors:
135: * ==================
136: *>
137: *> Ren-Cang Li, Computer Science Division, University of California
138: *> at Berkeley, USA
139: *>
140: * =====================================================================
141: SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
142: *
143: * -- LAPACK computational routine (version 3.6.0) --
144: * -- LAPACK is a software package provided by Univ. of Tennessee, --
145: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146: * November 2015
147: *
148: * .. Scalar Arguments ..
149: LOGICAL ORGATI
150: INTEGER INFO, KNITER
151: DOUBLE PRECISION FINIT, RHO, TAU
152: * ..
153: * .. Array Arguments ..
154: DOUBLE PRECISION D( 3 ), Z( 3 )
155: * ..
156: *
157: * =====================================================================
158: *
159: * .. Parameters ..
160: INTEGER MAXIT
161: PARAMETER ( MAXIT = 40 )
162: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT
163: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
164: $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 )
165: * ..
166: * .. External Functions ..
167: DOUBLE PRECISION DLAMCH
168: EXTERNAL DLAMCH
169: * ..
170: * .. Local Arrays ..
171: DOUBLE PRECISION DSCALE( 3 ), ZSCALE( 3 )
172: * ..
173: * .. Local Scalars ..
174: LOGICAL SCALE
175: INTEGER I, ITER, NITER
176: DOUBLE PRECISION A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
177: $ FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
178: $ SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4,
179: $ LBD, UBD
180: * ..
181: * .. Intrinsic Functions ..
182: INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT
183: * ..
184: * .. Executable Statements ..
185: *
186: INFO = 0
187: *
188: IF( ORGATI ) THEN
189: LBD = D(2)
190: UBD = D(3)
191: ELSE
192: LBD = D(1)
193: UBD = D(2)
194: END IF
195: IF( FINIT .LT. ZERO )THEN
196: LBD = ZERO
197: ELSE
198: UBD = ZERO
199: END IF
200: *
201: NITER = 1
202: TAU = ZERO
203: IF( KNITER.EQ.2 ) THEN
204: IF( ORGATI ) THEN
205: TEMP = ( D( 3 )-D( 2 ) ) / TWO
206: C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
207: A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
208: B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
209: ELSE
210: TEMP = ( D( 1 )-D( 2 ) ) / TWO
211: C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
212: A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
213: B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
214: END IF
215: TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
216: A = A / TEMP
217: B = B / TEMP
218: C = C / TEMP
219: IF( C.EQ.ZERO ) THEN
220: TAU = B / A
221: ELSE IF( A.LE.ZERO ) THEN
222: TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
223: ELSE
224: TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
225: END IF
226: IF( TAU .LT. LBD .OR. TAU .GT. UBD )
227: $ TAU = ( LBD+UBD )/TWO
228: IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN
229: TAU = ZERO
230: ELSE
231: TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) +
232: $ TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) +
233: $ TAU*Z(3)/( D(3)*( D( 3 )-TAU ) )
234: IF( TEMP .LE. ZERO )THEN
235: LBD = TAU
236: ELSE
237: UBD = TAU
238: END IF
239: IF( ABS( FINIT ).LE.ABS( TEMP ) )
240: $ TAU = ZERO
241: END IF
242: END IF
243: *
244: * get machine parameters for possible scaling to avoid overflow
245: *
246: * modified by Sven: parameters SMALL1, SMINV1, SMALL2,
247: * SMINV2, EPS are not SAVEd anymore between one call to the
248: * others but recomputed at each call
249: *
250: EPS = DLAMCH( 'Epsilon' )
251: BASE = DLAMCH( 'Base' )
252: SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) /
253: $ THREE ) )
254: SMINV1 = ONE / SMALL1
255: SMALL2 = SMALL1*SMALL1
256: SMINV2 = SMINV1*SMINV1
257: *
258: * Determine if scaling of inputs necessary to avoid overflow
259: * when computing 1/TEMP**3
260: *
261: IF( ORGATI ) THEN
262: TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
263: ELSE
264: TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
265: END IF
266: SCALE = .FALSE.
267: IF( TEMP.LE.SMALL1 ) THEN
268: SCALE = .TRUE.
269: IF( TEMP.LE.SMALL2 ) THEN
270: *
271: * Scale up by power of radix nearest 1/SAFMIN**(2/3)
272: *
273: SCLFAC = SMINV2
274: SCLINV = SMALL2
275: ELSE
276: *
277: * Scale up by power of radix nearest 1/SAFMIN**(1/3)
278: *
279: SCLFAC = SMINV1
280: SCLINV = SMALL1
281: END IF
282: *
283: * Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
284: *
285: DO 10 I = 1, 3
286: DSCALE( I ) = D( I )*SCLFAC
287: ZSCALE( I ) = Z( I )*SCLFAC
288: 10 CONTINUE
289: TAU = TAU*SCLFAC
290: LBD = LBD*SCLFAC
291: UBD = UBD*SCLFAC
292: ELSE
293: *
294: * Copy D and Z to DSCALE and ZSCALE
295: *
296: DO 20 I = 1, 3
297: DSCALE( I ) = D( I )
298: ZSCALE( I ) = Z( I )
299: 20 CONTINUE
300: END IF
301: *
302: FC = ZERO
303: DF = ZERO
304: DDF = ZERO
305: DO 30 I = 1, 3
306: TEMP = ONE / ( DSCALE( I )-TAU )
307: TEMP1 = ZSCALE( I )*TEMP
308: TEMP2 = TEMP1*TEMP
309: TEMP3 = TEMP2*TEMP
310: FC = FC + TEMP1 / DSCALE( I )
311: DF = DF + TEMP2
312: DDF = DDF + TEMP3
313: 30 CONTINUE
314: F = FINIT + TAU*FC
315: *
316: IF( ABS( F ).LE.ZERO )
317: $ GO TO 60
318: IF( F .LE. ZERO )THEN
319: LBD = TAU
320: ELSE
321: UBD = TAU
322: END IF
323: *
324: * Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
325: * scheme
326: *
327: * It is not hard to see that
328: *
329: * 1) Iterations will go up monotonically
330: * if FINIT < 0;
331: *
332: * 2) Iterations will go down monotonically
333: * if FINIT > 0.
334: *
335: ITER = NITER + 1
336: *
337: DO 50 NITER = ITER, MAXIT
338: *
339: IF( ORGATI ) THEN
340: TEMP1 = DSCALE( 2 ) - TAU
341: TEMP2 = DSCALE( 3 ) - TAU
342: ELSE
343: TEMP1 = DSCALE( 1 ) - TAU
344: TEMP2 = DSCALE( 2 ) - TAU
345: END IF
346: A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
347: B = TEMP1*TEMP2*F
348: C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
349: TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
350: A = A / TEMP
351: B = B / TEMP
352: C = C / TEMP
353: IF( C.EQ.ZERO ) THEN
354: ETA = B / A
355: ELSE IF( A.LE.ZERO ) THEN
356: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
357: ELSE
358: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
359: END IF
360: IF( F*ETA.GE.ZERO ) THEN
361: ETA = -F / DF
362: END IF
363: *
364: TAU = TAU + ETA
365: IF( TAU .LT. LBD .OR. TAU .GT. UBD )
366: $ TAU = ( LBD + UBD )/TWO
367: *
368: FC = ZERO
369: ERRETM = ZERO
370: DF = ZERO
371: DDF = ZERO
372: DO 40 I = 1, 3
373: IF ( ( DSCALE( I )-TAU ).NE.ZERO ) THEN
374: TEMP = ONE / ( DSCALE( I )-TAU )
375: TEMP1 = ZSCALE( I )*TEMP
376: TEMP2 = TEMP1*TEMP
377: TEMP3 = TEMP2*TEMP
378: TEMP4 = TEMP1 / DSCALE( I )
379: FC = FC + TEMP4
380: ERRETM = ERRETM + ABS( TEMP4 )
381: DF = DF + TEMP2
382: DDF = DDF + TEMP3
383: ELSE
384: GO TO 60
385: END IF
386: 40 CONTINUE
387: F = FINIT + TAU*FC
388: ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
389: $ ABS( TAU )*DF
390: IF( ( ABS( F ).LE.FOUR*EPS*ERRETM ) .OR.
391: $ ( (UBD-LBD).LE.FOUR*EPS*ABS(TAU) ) )
392: $ GO TO 60
393: IF( F .LE. ZERO )THEN
394: LBD = TAU
395: ELSE
396: UBD = TAU
397: END IF
398: 50 CONTINUE
399: INFO = 1
400: 60 CONTINUE
401: *
402: * Undo scaling
403: *
404: IF( SCALE )
405: $ TAU = TAU*SCLINV
406: RETURN
407: *
408: * End of DLAED6
409: *
410: 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