1*> \brief \b SLAED6 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 SLAED6 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed6.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed6.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed6.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
22*
23*       .. Scalar Arguments ..
24*       LOGICAL            ORGATI
25*       INTEGER            INFO, KNITER
26*       REAL               FINIT, RHO, TAU
27*       ..
28*       .. Array Arguments ..
29*       REAL               D( 3 ), Z( 3 )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> SLAED6 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 SLAED4 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 SLAED4 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*>               SLAED4 for further details.
69*> \endverbatim
70*>
71*> \param[in] RHO
72*> \verbatim
73*>          RHO is REAL
74*>               Refer to the equation f(x) above.
75*> \endverbatim
76*>
77*> \param[in] D
78*> \verbatim
79*>          D is REAL array, dimension (3)
80*>               D satisfies d(1) < d(2) < d(3).
81*> \endverbatim
82*>
83*> \param[in] Z
84*> \verbatim
85*>          Z is REAL 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 REAL
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 REAL
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 SLAED6( 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      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*
410      END
411