LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slaed6()

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.
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 slaed6.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  REAL FINIT, RHO, TAU
149 * ..
150 * .. Array Arguments ..
151  REAL D( 3 ), Z( 3 )
152 * ..
153 *
154 * =====================================================================
155 *
156 * .. Parameters ..
157  INTEGER MAXIT
158  parameter( maxit = 40 )
159  REAL ZERO, ONE, TWO, THREE, FOUR, EIGHT
160  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
161  $ three = 3.0e0, four = 4.0e0, eight = 8.0e0 )
162 * ..
163 * .. External Functions ..
164  REAL SLAMCH
165  EXTERNAL slamch
166 * ..
167 * .. Local Arrays ..
168  REAL DSCALE( 3 ), ZSCALE( 3 )
169 * ..
170 * .. Local Scalars ..
171  LOGICAL SCALE
172  INTEGER I, ITER, NITER
173  REAL 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 = slamch( 'Epsilon' )
248  base = slamch( 'Base' )
249  small1 = base**( int( log( slamch( '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 SLAED6
406 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the caller graph for this function: