LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine slaed6 ( integer  KNITER,
logical  ORGATI,
real  RHO,
real, dimension( 3 )  D,
real, dimension( 3 )  Z,
real  FINIT,
real  TAU,
integer  INFO 
)

SLAED6 used by sstedc. Computes one Newton step in solution of the secular equation.

Download SLAED6 + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 SLAED6 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 SLAED4 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 SLAED4 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
               SLAED4 for further details.
[in]RHO
          RHO is REAL
               Refer to the equation f(x) above.
[in]D
          D is REAL array, dimension (3)
               D satisfies d(1) < d(2) < d(3).
[in]Z
          Z is REAL array, dimension (3)
               Each of the elements in z must be positive.
[in]FINIT
          FINIT is REAL
               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 REAL
               The root of the equation f(x).
[out]INFO
          INFO is INTEGER
               = 0: successful exit
               > 0: if INFO = 1, failure to converge
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
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 142 of file slaed6.f.

142 *
143 * -- LAPACK computational routine (version 3.6.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 * November 2015
147 *
148 * .. Scalar Arguments ..
149  LOGICAL orgati
150  INTEGER info, kniter
151  REAL finit, rho, tau
152 * ..
153 * .. Array Arguments ..
154  REAL d( 3 ), z( 3 )
155 * ..
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160  INTEGER maxit
161  parameter ( maxit = 40 )
162  REAL zero, one, two, three, four, eight
163  parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
164  $ three = 3.0e0, four = 4.0e0, eight = 8.0e0 )
165 * ..
166 * .. External Functions ..
167  REAL slamch
168  EXTERNAL slamch
169 * ..
170 * .. Local Arrays ..
171  REAL dscale( 3 ), zscale( 3 )
172 * ..
173 * .. Local Scalars ..
174  LOGICAL scale
175  INTEGER i, iter, niter
176  REAL 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 = slamch( 'Epsilon' )
251  base = slamch( 'Base' )
252  small1 = base**( int( log( slamch( '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 SLAED6
409 *
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69

Here is the caller graph for this function: