1      SUBROUTINE SLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
2*
3*  -- LAPACK routine (instrumented to count operations, version 3.0) --
4*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
5*     Courant Institute, NAG Ltd., and Rice University
6*     June 30, 1999
7*
8*     .. Scalar Arguments ..
9      LOGICAL            ORGATI
10      INTEGER            INFO, KNITER
11      REAL               FINIT, RHO, TAU
12*     ..
13*     .. Array Arguments ..
14      REAL               D( 3 ), Z( 3 )
15*     ..
16*     Common block to return operation count and iteration count
17*     ITCNT is unchanged, OPS is only incremented
18*     .. Common blocks ..
19      COMMON             / LATIME / OPS, ITCNT
20*     ..
21*     .. Scalars in Common ..
22      REAL               ITCNT, OPS
23*     ..
24*
25*  Purpose
26*  =======
27*
28*  SLAED6 computes the positive or negative root (closest to the origin)
29*  of
30*                   z(1)        z(2)        z(3)
31*  f(x) =   rho + --------- + ---------- + ---------
32*                  d(1)-x      d(2)-x      d(3)-x
33*
34*  It is assumed that
35*
36*        if ORGATI = .true. the root is between d(2) and d(3);
37*        otherwise it is between d(1) and d(2)
38*
39*  This routine will be called by SLAED4 when necessary. In most cases,
40*  the root sought is the smallest in magnitude, though it might not be
41*  in some extremely rare situations.
42*
43*  Arguments
44*  =========
45*
46*  KNITER       (input) INTEGER
47*               Refer to SLAED4 for its significance.
48*
49*  ORGATI       (input) LOGICAL
50*               If ORGATI is true, the needed root is between d(2) and
51*               d(3); otherwise it is between d(1) and d(2).  See
52*               SLAED4 for further details.
53*
54*  RHO          (input) REAL
55*               Refer to the equation f(x) above.
56*
57*  D            (input) REAL array, dimension (3)
58*               D satisfies d(1) < d(2) < d(3).
59*
60*  Z            (input) REAL array, dimension (3)
61*               Each of the elements in z must be positive.
62*
63*  FINIT        (input) REAL
64*               The value of f at 0. It is more accurate than the one
65*               evaluated inside this routine (if someone wants to do
66*               so).
67*
68*  TAU          (output) REAL
69*               The root of the equation f(x).
70*
71*  INFO         (output) INTEGER
72*               = 0: successful exit
73*               > 0: if INFO = 1, failure to converge
74*
75*  Further Details
76*  ===============
77*
78*  Based on contributions by
79*     Ren-Cang Li, Computer Science Division, University of California
80*     at Berkeley, USA
81*
82*  =====================================================================
83*
84*     .. Parameters ..
85      INTEGER            MAXIT
86      PARAMETER          ( MAXIT = 20 )
87      REAL               ZERO, ONE, TWO, THREE, FOUR, EIGHT
88      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
89     $                   THREE = 3.0E0, FOUR = 4.0E0, EIGHT = 8.0E0 )
90*     ..
91*     .. External Functions ..
92      REAL               SLAMCH
93      EXTERNAL           SLAMCH
94*     ..
95*     .. Local Arrays ..
96      REAL               DSCALE( 3 ), ZSCALE( 3 )
97*     ..
98*     .. Local Scalars ..
99      LOGICAL            FIRST, SCALE
100      INTEGER            I, ITER, NITER
101      REAL               A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
102     $                   FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
103     $                   SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4
104*     ..
105*     .. Save statement ..
106      SAVE               FIRST, SMALL1, SMINV1, SMALL2, SMINV2, EPS
107*     ..
108*     .. Intrinsic Functions ..
109      INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
110*     ..
111*     .. Data statements ..
112      DATA               FIRST / .TRUE. /
113*     ..
114*     .. Executable Statements ..
115*
116      INFO = 0
117*
118      NITER = 1
119      TAU = ZERO
120      IF( KNITER.EQ.2 ) THEN
121         IF( ORGATI ) THEN
122            TEMP = ( D( 3 )-D( 2 ) ) / TWO
123            C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
124            A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
125            B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
126         ELSE
127            TEMP = ( D( 1 )-D( 2 ) ) / TWO
128            C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
129            A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
130            B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
131         END IF
132         TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
133         A = A / TEMP
134         B = B / TEMP
135         C = C / TEMP
136         OPS = OPS + 19
137         IF( C.EQ.ZERO ) THEN
138            TAU = B / A
139            OPS = OPS + 1
140         ELSE IF( A.LE.ZERO ) THEN
141            OPS = OPS + 8
142            TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
143         ELSE
144            OPS = OPS + 8
145            TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
146         END IF
147         OPS = OPS + 9
148         TEMP = RHO + Z( 1 ) / ( D( 1 )-TAU ) +
149     $          Z( 2 ) / ( D( 2 )-TAU ) + Z( 3 ) / ( D( 3 )-TAU )
150         IF( ABS( FINIT ).LE.ABS( TEMP ) )
151     $      TAU = ZERO
152      END IF
153*
154*     On first call to routine, get machine parameters for
155*     possible scaling to avoid overflow
156*
157      IF( FIRST ) THEN
158         EPS = SLAMCH( 'Epsilon' )
159         BASE = SLAMCH( 'Base' )
160         SMALL1 = BASE**( INT( LOG( SLAMCH( 'SafMin' ) ) / LOG( BASE ) /
161     $            THREE ) )
162         SMINV1 = ONE / SMALL1
163         SMALL2 = SMALL1*SMALL1
164         SMINV2 = SMINV1*SMINV1
165         FIRST = .FALSE.
166      END IF
167*
168*     Determine if scaling of inputs necessary to avoid overflow
169*     when computing 1/TEMP**3
170*
171      OPS = OPS + 2
172      IF( ORGATI ) THEN
173         TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
174      ELSE
175         TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
176      END IF
177      SCALE = .FALSE.
178      IF( TEMP.LE.SMALL1 ) THEN
179         SCALE = .TRUE.
180         IF( TEMP.LE.SMALL2 ) THEN
181*
182*        Scale up by power of radix nearest 1/SAFMIN**(2/3)
183*
184            SCLFAC = SMINV2
185            SCLINV = SMALL2
186         ELSE
187*
188*        Scale up by power of radix nearest 1/SAFMIN**(1/3)
189*
190            SCLFAC = SMINV1
191            SCLINV = SMALL1
192         END IF
193*
194*        Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
195*
196         OPS = OPS + 7
197         DO 10 I = 1, 3
198            DSCALE( I ) = D( I )*SCLFAC
199            ZSCALE( I ) = Z( I )*SCLFAC
200   10    CONTINUE
201         TAU = TAU*SCLFAC
202      ELSE
203*
204*        Copy D and Z to DSCALE and ZSCALE
205*
206         DO 20 I = 1, 3
207            DSCALE( I ) = D( I )
208            ZSCALE( I ) = Z( I )
209   20    CONTINUE
210      END IF
211*
212      FC = ZERO
213      DF = ZERO
214      DDF = ZERO
215      OPS = OPS + 11
216      DO 30 I = 1, 3
217         TEMP = ONE / ( DSCALE( I )-TAU )
218         TEMP1 = ZSCALE( I )*TEMP
219         TEMP2 = TEMP1*TEMP
220         TEMP3 = TEMP2*TEMP
221         FC = FC + TEMP1 / DSCALE( I )
222         DF = DF + TEMP2
223         DDF = DDF + TEMP3
224   30 CONTINUE
225      F = FINIT + TAU*FC
226*
227      IF( ABS( F ).LE.ZERO )
228     $   GO TO 60
229*
230*        Iteration begins
231*
232*     It is not hard to see that
233*
234*           1) Iterations will go up monotonically
235*              if FINIT < 0;
236*
237*           2) Iterations will go down monotonically
238*              if FINIT > 0.
239*
240      ITER = NITER + 1
241*
242      DO 50 NITER = ITER, MAXIT
243*
244         OPS = OPS + 18
245         IF( ORGATI ) THEN
246            TEMP1 = DSCALE( 2 ) - TAU
247            TEMP2 = DSCALE( 3 ) - TAU
248         ELSE
249            TEMP1 = DSCALE( 1 ) - TAU
250            TEMP2 = DSCALE( 2 ) - TAU
251         END IF
252         A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
253         B = TEMP1*TEMP2*F
254         C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
255         TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
256         A = A / TEMP
257         B = B / TEMP
258         C = C / TEMP
259         IF( C.EQ.ZERO ) THEN
260            OPS = OPS + 1
261            ETA = B / A
262         ELSE IF( A.LE.ZERO ) THEN
263            OPS = OPS + 8
264            ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
265         ELSE
266            OPS = OPS + 8
267            ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
268         END IF
269         IF( F*ETA.GE.ZERO ) THEN
270            OPS = OPS + 1
271            ETA = -F / DF
272         END IF
273*
274         OPS = OPS + 1
275         TEMP = ETA + TAU
276         IF( ORGATI ) THEN
277            IF( ETA.GT.ZERO .AND. TEMP.GE.DSCALE( 3 ) ) THEN
278               OPS = OPS + 2
279               ETA = ( DSCALE( 3 )-TAU ) / TWO
280            END IF
281            IF( ETA.LT.ZERO .AND. TEMP.LE.DSCALE( 2 ) ) THEN
282               OPS = OPS + 2
283               ETA = ( DSCALE( 2 )-TAU ) / TWO
284            END IF
285         ELSE
286            IF( ETA.GT.ZERO .AND. TEMP.GE.DSCALE( 2 ) ) THEN
287               OPS = OPS + 2
288               ETA = ( DSCALE( 2 )-TAU ) / TWO
289            END IF
290            IF( ETA.LT.ZERO .AND. TEMP.LE.DSCALE( 1 ) ) THEN
291               OPS = OPS + 2
292               ETA = ( DSCALE( 1 )-TAU ) / TWO
293            END IF
294         END IF
295         OPS = OPS + 1
296         TAU = TAU + ETA
297*
298         FC = ZERO
299         ERRETM = ZERO
300         DF = ZERO
301         DDF = ZERO
302         OPS = OPS + 37
303         DO 40 I = 1, 3
304            TEMP = ONE / ( DSCALE( I )-TAU )
305            TEMP1 = ZSCALE( I )*TEMP
306            TEMP2 = TEMP1*TEMP
307            TEMP3 = TEMP2*TEMP
308            TEMP4 = TEMP1 / DSCALE( I )
309            FC = FC + TEMP4
310            ERRETM = ERRETM + ABS( TEMP4 )
311            DF = DF + TEMP2
312            DDF = DDF + TEMP3
313   40    CONTINUE
314         F = FINIT + TAU*FC
315         ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
316     $            ABS( TAU )*DF
317         IF( ABS( F ).LE.EPS*ERRETM )
318     $      GO TO 60
319   50 CONTINUE
320      INFO = 1
321   60 CONTINUE
322*
323*     Undo scaling
324*
325      IF( SCALE ) THEN
326         OPS = OPS + 1
327         TAU = TAU*SCLINV
328      END IF
329      RETURN
330*
331*     End of SLAED6
332*
333      END
334