LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaed6.f
Go to the documentation of this file.
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 laed6
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
subroutine dlaed6(kniter, orgati, rho, d, z, finit, tau, info)
DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation.
Definition dlaed6.f:140