LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dlaed6()

 subroutine dlaed6 ( integer KNITER, logical ORGATI, double precision RHO, double precision, dimension( 3 ) D, double precision, dimension( 3 ) Z, double precision FINIT, double precision TAU, integer INFO )

DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation.

Purpose:
DLAED6 computes the positive or negative root (closest to the origin)
of
z(1)        z(2)        z(3)
f(x) =   rho + --------- + ---------- + ---------
d(1)-x      d(2)-x      d(3)-x

It is assumed that

if ORGATI = .true. the root is between d(2) and d(3);
otherwise it is between d(1) and d(2)

This routine will be called by DLAED4 when necessary. In most cases,
the root sought is the smallest in magnitude, though it might not be
in some extremely rare situations.
Parameters
 [in] KNITER KNITER is INTEGER Refer to DLAED4 for its significance. [in] ORGATI ORGATI is LOGICAL If ORGATI is true, the needed root is between d(2) and d(3); otherwise it is between d(1) and d(2). See DLAED4 for further details. [in] RHO RHO is DOUBLE PRECISION Refer to the equation f(x) above. [in] D D is DOUBLE PRECISION array, dimension (3) D satisfies d(1) < d(2) < d(3). [in] Z Z is DOUBLE PRECISION array, dimension (3) Each of the elements in z must be positive. [in] FINIT FINIT is DOUBLE PRECISION The value of f at 0. It is more accurate than the one evaluated inside this routine (if someone wants to do so). [out] TAU TAU is DOUBLE PRECISION The root of the equation f(x). [out] INFO INFO is INTEGER = 0: successful exit > 0: if INFO = 1, failure to converge
Further Details:
10/02/03: This version has a few statements commented out for thread
safety (machine parameters are computed on each entry). SJH.

05/10/06: Modified from a new version of Ren-Cang Li, use
Gragg-Thornton-Warner cubic convergent scheme for better stability.
Contributors:
Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

Definition at line 139 of file dlaed6.f.

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*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
Here is the caller graph for this function: