Annotation of rpl/lapack/lapack/dlaed6.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLAED6
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLAED6 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed6.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed6.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed6.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * LOGICAL ORGATI
! 25: * INTEGER INFO, KNITER
! 26: * DOUBLE PRECISION FINIT, RHO, TAU
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * DOUBLE PRECISION D( 3 ), Z( 3 )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> DLAED6 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 DLAED4 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 DLAED4 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: *> DLAED4 for further details.
! 69: *> \endverbatim
! 70: *>
! 71: *> \param[in] RHO
! 72: *> \verbatim
! 73: *> RHO is DOUBLE PRECISION
! 74: *> Refer to the equation f(x) above.
! 75: *> \endverbatim
! 76: *>
! 77: *> \param[in] D
! 78: *> \verbatim
! 79: *> D is DOUBLE PRECISION array, dimension (3)
! 80: *> D satisfies d(1) < d(2) < d(3).
! 81: *> \endverbatim
! 82: *>
! 83: *> \param[in] Z
! 84: *> \verbatim
! 85: *> Z is DOUBLE PRECISION 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 DOUBLE PRECISION
! 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 DOUBLE PRECISION
! 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 November 2011
! 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: * =====================================================================
1.1 bertrand 141: SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
142: *
1.9 ! bertrand 143: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 144: * -- LAPACK is a software package provided by Univ. of Tennessee, --
145: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 146: * November 2011
1.1 bertrand 147: *
148: * .. Scalar Arguments ..
149: LOGICAL ORGATI
150: INTEGER INFO, KNITER
151: DOUBLE PRECISION FINIT, RHO, TAU
152: * ..
153: * .. Array Arguments ..
154: DOUBLE PRECISION D( 3 ), Z( 3 )
155: * ..
156: *
157: * =====================================================================
158: *
159: * .. Parameters ..
160: INTEGER MAXIT
161: PARAMETER ( MAXIT = 40 )
162: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT
163: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
164: $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 )
165: * ..
166: * .. External Functions ..
167: DOUBLE PRECISION DLAMCH
168: EXTERNAL DLAMCH
169: * ..
170: * .. Local Arrays ..
171: DOUBLE PRECISION DSCALE( 3 ), ZSCALE( 3 )
172: * ..
173: * .. Local Scalars ..
174: LOGICAL SCALE
175: INTEGER I, ITER, NITER
176: DOUBLE PRECISION 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 = DLAMCH( 'Epsilon' )
251: BASE = DLAMCH( 'Base' )
252: SMALL1 = BASE**( INT( LOG( DLAMCH( '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: TEMP = ONE / ( DSCALE( I )-TAU )
374: TEMP1 = ZSCALE( I )*TEMP
375: TEMP2 = TEMP1*TEMP
376: TEMP3 = TEMP2*TEMP
377: TEMP4 = TEMP1 / DSCALE( I )
378: FC = FC + TEMP4
379: ERRETM = ERRETM + ABS( TEMP4 )
380: DF = DF + TEMP2
381: DDF = DDF + TEMP3
382: 40 CONTINUE
383: F = FINIT + TAU*FC
384: ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
385: $ ABS( TAU )*DF
386: IF( ABS( F ).LE.EPS*ERRETM )
387: $ GO TO 60
388: IF( F .LE. ZERO )THEN
389: LBD = TAU
390: ELSE
391: UBD = TAU
392: END IF
393: 50 CONTINUE
394: INFO = 1
395: 60 CONTINUE
396: *
397: * Undo scaling
398: *
399: IF( SCALE )
400: $ TAU = TAU*SCLINV
401: RETURN
402: *
403: * End of DLAED6
404: *
405: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed6.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack