![[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|
Return to dlaed6.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * February 2007
7: *
8: * .. Scalar Arguments ..
9: LOGICAL ORGATI
10: INTEGER INFO, KNITER
11: DOUBLE PRECISION FINIT, RHO, TAU
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION D( 3 ), Z( 3 )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DLAED6 computes the positive or negative root (closest to the origin)
21: * of
22: * z(1) z(2) z(3)
23: * f(x) = rho + --------- + ---------- + ---------
24: * d(1)-x d(2)-x d(3)-x
25: *
26: * It is assumed that
27: *
28: * if ORGATI = .true. the root is between d(2) and d(3);
29: * otherwise it is between d(1) and d(2)
30: *
31: * This routine will be called by DLAED4 when necessary. In most cases,
32: * the root sought is the smallest in magnitude, though it might not be
33: * in some extremely rare situations.
34: *
35: * Arguments
36: * =========
37: *
38: * KNITER (input) INTEGER
39: * Refer to DLAED4 for its significance.
40: *
41: * ORGATI (input) LOGICAL
42: * If ORGATI is true, the needed root is between d(2) and
43: * d(3); otherwise it is between d(1) and d(2). See
44: * DLAED4 for further details.
45: *
46: * RHO (input) DOUBLE PRECISION
47: * Refer to the equation f(x) above.
48: *
49: * D (input) DOUBLE PRECISION array, dimension (3)
50: * D satisfies d(1) < d(2) < d(3).
51: *
52: * Z (input) DOUBLE PRECISION array, dimension (3)
53: * Each of the elements in z must be positive.
54: *
55: * FINIT (input) DOUBLE PRECISION
56: * The value of f at 0. It is more accurate than the one
57: * evaluated inside this routine (if someone wants to do
58: * so).
59: *
60: * TAU (output) DOUBLE PRECISION
61: * The root of the equation f(x).
62: *
63: * INFO (output) INTEGER
64: * = 0: successful exit
65: * > 0: if INFO = 1, failure to converge
66: *
67: * Further Details
68: * ===============
69: *
70: * 30/06/99: Based on contributions by
71: * Ren-Cang Li, Computer Science Division, University of California
72: * at Berkeley, USA
73: *
74: * 10/02/03: This version has a few statements commented out for thread
75: * safety (machine parameters are computed on each entry). SJH.
76: *
77: * 05/10/06: Modified from a new version of Ren-Cang Li, use
78: * Gragg-Thornton-Warner cubic convergent scheme for better stability.
79: *
80: * =====================================================================
81: *
82: * .. Parameters ..
83: INTEGER MAXIT
84: PARAMETER ( MAXIT = 40 )
85: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT
86: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
87: $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 )
88: * ..
89: * .. External Functions ..
90: DOUBLE PRECISION DLAMCH
91: EXTERNAL DLAMCH
92: * ..
93: * .. Local Arrays ..
94: DOUBLE PRECISION DSCALE( 3 ), ZSCALE( 3 )
95: * ..
96: * .. Local Scalars ..
97: LOGICAL SCALE
98: INTEGER I, ITER, NITER
99: DOUBLE PRECISION A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
100: $ FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
101: $ SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4,
102: $ LBD, UBD
103: * ..
104: * .. Intrinsic Functions ..
105: INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT
106: * ..
107: * .. Executable Statements ..
108: *
109: INFO = 0
110: *
111: IF( ORGATI ) THEN
112: LBD = D(2)
113: UBD = D(3)
114: ELSE
115: LBD = D(1)
116: UBD = D(2)
117: END IF
118: IF( FINIT .LT. ZERO )THEN
119: LBD = ZERO
120: ELSE
121: UBD = ZERO
122: END IF
123: *
124: NITER = 1
125: TAU = ZERO
126: IF( KNITER.EQ.2 ) THEN
127: IF( ORGATI ) THEN
128: TEMP = ( D( 3 )-D( 2 ) ) / TWO
129: C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
130: A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
131: B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
132: ELSE
133: TEMP = ( D( 1 )-D( 2 ) ) / TWO
134: C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
135: A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
136: B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
137: END IF
138: TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
139: A = A / TEMP
140: B = B / TEMP
141: C = C / TEMP
142: IF( C.EQ.ZERO ) THEN
143: TAU = B / A
144: ELSE IF( A.LE.ZERO ) THEN
145: TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
146: ELSE
147: TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
148: END IF
149: IF( TAU .LT. LBD .OR. TAU .GT. UBD )
150: $ TAU = ( LBD+UBD )/TWO
151: IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN
152: TAU = ZERO
153: ELSE
154: TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) +
155: $ TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) +
156: $ TAU*Z(3)/( D(3)*( D( 3 )-TAU ) )
157: IF( TEMP .LE. ZERO )THEN
158: LBD = TAU
159: ELSE
160: UBD = TAU
161: END IF
162: IF( ABS( FINIT ).LE.ABS( TEMP ) )
163: $ TAU = ZERO
164: END IF
165: END IF
166: *
167: * get machine parameters for possible scaling to avoid overflow
168: *
169: * modified by Sven: parameters SMALL1, SMINV1, SMALL2,
170: * SMINV2, EPS are not SAVEd anymore between one call to the
171: * others but recomputed at each call
172: *
173: EPS = DLAMCH( 'Epsilon' )
174: BASE = DLAMCH( 'Base' )
175: SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) /
176: $ THREE ) )
177: SMINV1 = ONE / SMALL1
178: SMALL2 = SMALL1*SMALL1
179: SMINV2 = SMINV1*SMINV1
180: *
181: * Determine if scaling of inputs necessary to avoid overflow
182: * when computing 1/TEMP**3
183: *
184: IF( ORGATI ) THEN
185: TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
186: ELSE
187: TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
188: END IF
189: SCALE = .FALSE.
190: IF( TEMP.LE.SMALL1 ) THEN
191: SCALE = .TRUE.
192: IF( TEMP.LE.SMALL2 ) THEN
193: *
194: * Scale up by power of radix nearest 1/SAFMIN**(2/3)
195: *
196: SCLFAC = SMINV2
197: SCLINV = SMALL2
198: ELSE
199: *
200: * Scale up by power of radix nearest 1/SAFMIN**(1/3)
201: *
202: SCLFAC = SMINV1
203: SCLINV = SMALL1
204: END IF
205: *
206: * Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
207: *
208: DO 10 I = 1, 3
209: DSCALE( I ) = D( I )*SCLFAC
210: ZSCALE( I ) = Z( I )*SCLFAC
211: 10 CONTINUE
212: TAU = TAU*SCLFAC
213: LBD = LBD*SCLFAC
214: UBD = UBD*SCLFAC
215: ELSE
216: *
217: * Copy D and Z to DSCALE and ZSCALE
218: *
219: DO 20 I = 1, 3
220: DSCALE( I ) = D( I )
221: ZSCALE( I ) = Z( I )
222: 20 CONTINUE
223: END IF
224: *
225: FC = ZERO
226: DF = ZERO
227: DDF = ZERO
228: DO 30 I = 1, 3
229: TEMP = ONE / ( DSCALE( I )-TAU )
230: TEMP1 = ZSCALE( I )*TEMP
231: TEMP2 = TEMP1*TEMP
232: TEMP3 = TEMP2*TEMP
233: FC = FC + TEMP1 / DSCALE( I )
234: DF = DF + TEMP2
235: DDF = DDF + TEMP3
236: 30 CONTINUE
237: F = FINIT + TAU*FC
238: *
239: IF( ABS( F ).LE.ZERO )
240: $ GO TO 60
241: IF( F .LE. ZERO )THEN
242: LBD = TAU
243: ELSE
244: UBD = TAU
245: END IF
246: *
247: * Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
248: * scheme
249: *
250: * It is not hard to see that
251: *
252: * 1) Iterations will go up monotonically
253: * if FINIT < 0;
254: *
255: * 2) Iterations will go down monotonically
256: * if FINIT > 0.
257: *
258: ITER = NITER + 1
259: *
260: DO 50 NITER = ITER, MAXIT
261: *
262: IF( ORGATI ) THEN
263: TEMP1 = DSCALE( 2 ) - TAU
264: TEMP2 = DSCALE( 3 ) - TAU
265: ELSE
266: TEMP1 = DSCALE( 1 ) - TAU
267: TEMP2 = DSCALE( 2 ) - TAU
268: END IF
269: A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
270: B = TEMP1*TEMP2*F
271: C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
272: TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
273: A = A / TEMP
274: B = B / TEMP
275: C = C / TEMP
276: IF( C.EQ.ZERO ) THEN
277: ETA = B / A
278: ELSE IF( A.LE.ZERO ) THEN
279: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
280: ELSE
281: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
282: END IF
283: IF( F*ETA.GE.ZERO ) THEN
284: ETA = -F / DF
285: END IF
286: *
287: TAU = TAU + ETA
288: IF( TAU .LT. LBD .OR. TAU .GT. UBD )
289: $ TAU = ( LBD + UBD )/TWO
290: *
291: FC = ZERO
292: ERRETM = ZERO
293: DF = ZERO
294: DDF = ZERO
295: DO 40 I = 1, 3
296: TEMP = ONE / ( DSCALE( I )-TAU )
297: TEMP1 = ZSCALE( I )*TEMP
298: TEMP2 = TEMP1*TEMP
299: TEMP3 = TEMP2*TEMP
300: TEMP4 = TEMP1 / DSCALE( I )
301: FC = FC + TEMP4
302: ERRETM = ERRETM + ABS( TEMP4 )
303: DF = DF + TEMP2
304: DDF = DDF + TEMP3
305: 40 CONTINUE
306: F = FINIT + TAU*FC
307: ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
308: $ ABS( TAU )*DF
309: IF( ABS( F ).LE.EPS*ERRETM )
310: $ GO TO 60
311: IF( F .LE. ZERO )THEN
312: LBD = TAU
313: ELSE
314: UBD = TAU
315: END IF
316: 50 CONTINUE
317: INFO = 1
318: 60 CONTINUE
319: *
320: * Undo scaling
321: *
322: IF( SCALE )
323: $ TAU = TAU*SCLINV
324: RETURN
325: *
326: * End of DLAED6
327: *
328: END