001:       SUBROUTINE SLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     February 2007
006: *
007: *     .. Scalar Arguments ..
008:       LOGICAL            ORGATI
009:       INTEGER            INFO, KNITER
010:       REAL               FINIT, RHO, TAU
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               D( 3 ), Z( 3 )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  SLAED6 computes the positive or negative root (closest to the origin)
020: *  of
021: *                   z(1)        z(2)        z(3)
022: *  f(x) =   rho + --------- + ---------- + ---------
023: *                  d(1)-x      d(2)-x      d(3)-x
024: *
025: *  It is assumed that
026: *
027: *        if ORGATI = .true. the root is between d(2) and d(3);
028: *        otherwise it is between d(1) and d(2)
029: *
030: *  This routine will be called by SLAED4 when necessary. In most cases,
031: *  the root sought is the smallest in magnitude, though it might not be
032: *  in some extremely rare situations.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  KNITER       (input) INTEGER
038: *               Refer to SLAED4 for its significance.
039: *
040: *  ORGATI       (input) LOGICAL
041: *               If ORGATI is true, the needed root is between d(2) and
042: *               d(3); otherwise it is between d(1) and d(2).  See
043: *               SLAED4 for further details.
044: *
045: *  RHO          (input) REAL            
046: *               Refer to the equation f(x) above.
047: *
048: *  D            (input) REAL array, dimension (3)
049: *               D satisfies d(1) < d(2) < d(3).
050: *
051: *  Z            (input) REAL array, dimension (3)
052: *               Each of the elements in z must be positive.
053: *
054: *  FINIT        (input) REAL            
055: *               The value of f at 0. It is more accurate than the one
056: *               evaluated inside this routine (if someone wants to do
057: *               so).
058: *
059: *  TAU          (output) REAL            
060: *               The root of the equation f(x).
061: *
062: *  INFO         (output) INTEGER
063: *               = 0: successful exit
064: *               > 0: if INFO = 1, failure to converge
065: *
066: *  Further Details
067: *  ===============
068: *
069: *  30/06/99: Based on contributions by
070: *     Ren-Cang Li, Computer Science Division, University of California
071: *     at Berkeley, USA
072: *
073: *  10/02/03: This version has a few statements commented out for thread safety
074: *     (machine parameters are computed on each entry). SJH.
075: *
076: *  05/10/06: Modified from a new version of Ren-Cang Li, use
077: *     Gragg-Thornton-Warner cubic convergent scheme for better stability.
078: *
079: *  =====================================================================
080: *
081: *     .. Parameters ..
082:       INTEGER            MAXIT
083:       PARAMETER          ( MAXIT = 40 )
084:       REAL               ZERO, ONE, TWO, THREE, FOUR, EIGHT
085:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
086:      $                   THREE = 3.0E0, FOUR = 4.0E0, EIGHT = 8.0E0 )
087: *     ..
088: *     .. External Functions ..
089:       REAL               SLAMCH
090:       EXTERNAL           SLAMCH
091: *     ..
092: *     .. Local Arrays ..
093:       REAL               DSCALE( 3 ), ZSCALE( 3 )
094: *     ..
095: *     .. Local Scalars ..
096:       LOGICAL            SCALE
097:       INTEGER            I, ITER, NITER
098:       REAL               A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
099:      $                   FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
100:      $                   SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4, 
101:      $                   LBD, UBD
102: *     ..
103: *     .. Intrinsic Functions ..
104:       INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
105: *     ..
106: *     .. Executable Statements ..
107: *
108:       INFO = 0
109: *
110:       IF( ORGATI ) THEN
111:          LBD = D(2)
112:          UBD = D(3)
113:       ELSE
114:          LBD = D(1)
115:          UBD = D(2)
116:       END IF
117:       IF( FINIT .LT. ZERO )THEN
118:          LBD = ZERO
119:       ELSE
120:          UBD = ZERO 
121:       END IF
122: *
123:       NITER = 1
124:       TAU = ZERO
125:       IF( KNITER.EQ.2 ) THEN
126:          IF( ORGATI ) THEN
127:             TEMP = ( D( 3 )-D( 2 ) ) / TWO
128:             C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
129:             A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
130:             B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
131:          ELSE
132:             TEMP = ( D( 1 )-D( 2 ) ) / TWO
133:             C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
134:             A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
135:             B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
136:          END IF
137:          TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
138:          A = A / TEMP
139:          B = B / TEMP
140:          C = C / TEMP
141:          IF( C.EQ.ZERO ) THEN
142:             TAU = B / A
143:          ELSE IF( A.LE.ZERO ) THEN
144:             TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
145:          ELSE
146:             TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
147:          END IF
148:          IF( TAU .LT. LBD .OR. TAU .GT. UBD )
149:      $      TAU = ( LBD+UBD )/TWO
150:          IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN
151:             TAU = ZERO
152:          ELSE
153:             TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) +
154:      $                     TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) +
155:      $                     TAU*Z(3)/( D(3)*( D( 3 )-TAU ) )
156:             IF( TEMP .LE. ZERO )THEN
157:                LBD = TAU
158:             ELSE
159:                UBD = TAU
160:             END IF
161:             IF( ABS( FINIT ).LE.ABS( TEMP ) )
162:      $         TAU = ZERO
163:          END IF
164:       END IF
165: *
166: *     get machine parameters for possible scaling to avoid overflow
167: *
168: *     modified by Sven: parameters SMALL1, SMINV1, SMALL2,
169: *     SMINV2, EPS are not SAVEd anymore between one call to the
170: *     others but recomputed at each call
171: *
172:       EPS = SLAMCH( 'Epsilon' )
173:       BASE = SLAMCH( 'Base' )
174:       SMALL1 = BASE**( INT( LOG( SLAMCH( 'SafMin' ) ) / LOG( BASE ) /
175:      $         THREE ) )
176:       SMINV1 = ONE / SMALL1
177:       SMALL2 = SMALL1*SMALL1
178:       SMINV2 = SMINV1*SMINV1
179: *
180: *     Determine if scaling of inputs necessary to avoid overflow
181: *     when computing 1/TEMP**3
182: *
183:       IF( ORGATI ) THEN
184:          TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
185:       ELSE
186:          TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
187:       END IF
188:       SCALE = .FALSE.
189:       IF( TEMP.LE.SMALL1 ) THEN
190:          SCALE = .TRUE.
191:          IF( TEMP.LE.SMALL2 ) THEN
192: *
193: *        Scale up by power of radix nearest 1/SAFMIN**(2/3)
194: *
195:             SCLFAC = SMINV2
196:             SCLINV = SMALL2
197:          ELSE
198: *
199: *        Scale up by power of radix nearest 1/SAFMIN**(1/3)
200: *
201:             SCLFAC = SMINV1
202:             SCLINV = SMALL1
203:          END IF
204: *
205: *        Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
206: *
207:          DO 10 I = 1, 3
208:             DSCALE( I ) = D( I )*SCLFAC
209:             ZSCALE( I ) = Z( I )*SCLFAC
210:    10    CONTINUE
211:          TAU = TAU*SCLFAC
212:          LBD = LBD*SCLFAC
213:          UBD = UBD*SCLFAC
214:       ELSE
215: *
216: *        Copy D and Z to DSCALE and ZSCALE
217: *
218:          DO 20 I = 1, 3
219:             DSCALE( I ) = D( I )
220:             ZSCALE( I ) = Z( I )
221:    20    CONTINUE
222:       END IF
223: *
224:       FC = ZERO
225:       DF = ZERO
226:       DDF = ZERO
227:       DO 30 I = 1, 3
228:          TEMP = ONE / ( DSCALE( I )-TAU )
229:          TEMP1 = ZSCALE( I )*TEMP
230:          TEMP2 = TEMP1*TEMP
231:          TEMP3 = TEMP2*TEMP
232:          FC = FC + TEMP1 / DSCALE( I )
233:          DF = DF + TEMP2
234:          DDF = DDF + TEMP3
235:    30 CONTINUE
236:       F = FINIT + TAU*FC
237: *
238:       IF( ABS( F ).LE.ZERO )
239:      $   GO TO 60
240:       IF( F .LE. ZERO )THEN
241:          LBD = TAU
242:       ELSE
243:          UBD = TAU
244:       END IF
245: *
246: *        Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
247: *                            scheme
248: *
249: *     It is not hard to see that
250: *
251: *           1) Iterations will go up monotonically
252: *              if FINIT < 0;
253: *
254: *           2) Iterations will go down monotonically
255: *              if FINIT > 0.
256: *
257:       ITER = NITER + 1
258: *
259:       DO 50 NITER = ITER, MAXIT
260: *
261:          IF( ORGATI ) THEN
262:             TEMP1 = DSCALE( 2 ) - TAU
263:             TEMP2 = DSCALE( 3 ) - TAU
264:          ELSE
265:             TEMP1 = DSCALE( 1 ) - TAU
266:             TEMP2 = DSCALE( 2 ) - TAU
267:          END IF
268:          A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
269:          B = TEMP1*TEMP2*F
270:          C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
271:          TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
272:          A = A / TEMP
273:          B = B / TEMP
274:          C = C / TEMP
275:          IF( C.EQ.ZERO ) THEN
276:             ETA = B / A
277:          ELSE IF( A.LE.ZERO ) THEN
278:             ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
279:          ELSE
280:             ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
281:          END IF
282:          IF( F*ETA.GE.ZERO ) THEN
283:             ETA = -F / DF
284:          END IF
285: *
286:          TAU = TAU + ETA
287:          IF( TAU .LT. LBD .OR. TAU .GT. UBD )
288:      $      TAU = ( LBD + UBD )/TWO 
289: *
290:          FC = ZERO
291:          ERRETM = ZERO
292:          DF = ZERO
293:          DDF = ZERO
294:          DO 40 I = 1, 3
295:             TEMP = ONE / ( DSCALE( I )-TAU )
296:             TEMP1 = ZSCALE( I )*TEMP
297:             TEMP2 = TEMP1*TEMP
298:             TEMP3 = TEMP2*TEMP
299:             TEMP4 = TEMP1 / DSCALE( I )
300:             FC = FC + TEMP4
301:             ERRETM = ERRETM + ABS( TEMP4 )
302:             DF = DF + TEMP2
303:             DDF = DDF + TEMP3
304:    40    CONTINUE
305:          F = FINIT + TAU*FC
306:          ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
307:      $            ABS( TAU )*DF
308:          IF( ABS( F ).LE.EPS*ERRETM )
309:      $      GO TO 60
310:          IF( F .LE. ZERO )THEN
311:             LBD = TAU
312:          ELSE
313:             UBD = TAU
314:          END IF
315:    50 CONTINUE
316:       INFO = 1
317:    60 CONTINUE
318: *
319: *     Undo scaling
320: *
321:       IF( SCALE )
322:      $   TAU = TAU*SCLINV
323:       RETURN
324: *
325: *     End of SLAED6
326: *
327:       END
328: