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