LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 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
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