LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
dlaed6.f
Go to the documentation of this file.
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 December 2016
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.7.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 * December 2016
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
subroutine dlaed6(KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)
DLAED6 used by sstedc. Computes one Newton step in solution of the secular equation.
Definition: dlaed6.f:142