From 2f64e5a9cd832b9f837b9d48be9d51300bca5e49 Mon Sep 17 00:00:00 2001 From: julie Date: Tue, 11 Dec 2012 19:01:44 +0000 Subject: [PATCH] Commit Victor Liu's suggestion about making the complex division routine more robust tests and builds seems fine. Message sent on Oct 18th I just saw a paper on ArXiv about making the complex division routine more robust: http://arxiv.org/abs/1210.4539 Second author is actually the original inventor of the current algorithm in Lapack. I have attached my modified DLADIV routine, which passes all the tests in the build process. --- SRC/dladiv.f | 157 +++++++++++++++++++++++++++++++++++++++++++++++++++++------ SRC/sladiv.f | 157 +++++++++++++++++++++++++++++++++++++++++++++++++++++------ 2 files changed, 282 insertions(+), 32 deletions(-) diff --git a/SRC/dladiv.f b/SRC/dladiv.f index 306a6b0..c22d56d 100644 --- a/SRC/dladiv.f +++ b/SRC/dladiv.f @@ -36,8 +36,9 @@ *> p + i*q = --------- *> c + i*d *> -*> The algorithm is due to Robert L. Smith and can be found -*> in D. Knuth, The art of Computer Programming, Vol.2, p.195 +*> The algorithm is due to Michael Baudin and Robert L. Smith +*> and can be found in the paper +*> "A Robust Complex Division in Scilab" *> \endverbatim * * Arguments: @@ -83,17 +84,17 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date September 2012 +*> \date January 2013 * *> \ingroup auxOTHERauxiliary * * ===================================================================== SUBROUTINE DLADIV( A, B, C, D, P, Q ) * -* -- LAPACK auxiliary routine (version 3.4.2) -- +* -- LAPACK auxiliary routine (version 3.5.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* September 2012 +* January 2013 * * .. Scalar Arguments .. DOUBLE PRECISION A, B, C, D, P, Q @@ -101,28 +102,152 @@ * * ===================================================================== * +* .. Parameters .. + DOUBLE PRECISION BS + PARAMETER ( BS = 2.0D0 ) + DOUBLE PRECISION HALF + PARAMETER ( HALF = 0.5D0 ) + DOUBLE PRECISION TWO + PARAMETER ( TWO = 2.0D0 ) +* * .. Local Scalars .. - DOUBLE PRECISION E, F + DOUBLE PRECISION AA, BB, CC, DD, AB, CD, S, OV, UN, BE, EPS +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH + EXTERNAL DLAMCH +* .. +* .. External Subroutines .. + EXTERNAL DLADIV1 * .. * .. Intrinsic Functions .. - INTRINSIC ABS + INTRINSIC ABS, MAX * .. * .. Executable Statements .. * - IF( ABS( D ).LT.ABS( C ) ) THEN - E = D / C - F = C + D*E - P = ( A+B*E ) / F - Q = ( B-A*E ) / F + AA = A + BB = B + CC = C + DD = D + AB = MAX( ABS(A), ABS(B) ) + CD = MAX( ABS(C), ABS(D) ) + S = 1.0D0 + + OV = DLAMCH( 'Overflow threshold' ) + UN = DLAMCH( 'Safe minimum' ) + EPS = DLAMCH( 'Epsilon' ) + BE = BS / (EPS*EPS) + + IF( AB >= HALF*OV ) THEN + AA = HALF * AA + BB = HALF * BB + S = TWO * S + END IF + IF( CD >= HALF*OV ) THEN + CC = HALF * CC + DD = HALF * DD + S = HALF * S + END IF + IF( AB <= UN*BS/EPS ) THEN + AA = AA * BE + BB = BB * BE + S = S / BE + END IF + IF( CD <= UN*BS/EPS ) THEN + CC = CC * BE + DD = DD * BE + S = S * BE + END IF + IF( ABS( D ).LE.ABS( C ) ) THEN + CALL DLADIV1(AA, BB, CC, DD, P, Q) ELSE - E = C / D - F = D + C*E - P = ( B+A*E ) / F - Q = ( -A+B*E ) / F + CALL DLADIV1(BB, AA, DD, CC, P, Q) + Q = -Q END IF + P = P * S + Q = Q * S * RETURN * * End of DLADIV * END + + + + SUBROUTINE DLADIV1( A, B, C, D, P, Q ) +* +* -- LAPACK auxiliary routine (version 3.5.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* January 2013 +* +* .. Scalar Arguments .. + DOUBLE PRECISION A, B, C, D, P, Q +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D0 ) +* +* .. Local Scalars .. + DOUBLE PRECISION R, T +* .. +* .. External Functions .. + DOUBLE PRECISION DLADIV2 + EXTERNAL DLADIV2 +* .. +* .. Executable Statements .. +* + R = D / C + T = ONE / (C + D * R) + P = DLADIV2(A, B, C, D, R, T) + A = -A + Q = DLADIV2(B, A, C, D, R, T) +* + RETURN +* +* End of DLADIV1 +* + END + + DOUBLE PRECISION FUNCTION DLADIV2( A, B, C, D, R, T ) +* +* -- LAPACK auxiliary routine (version 3.5.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* January 2013 +* +* .. Scalar Arguments .. + DOUBLE PRECISION A, B, C, D, R, T +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D0 ) +* +* .. Local Scalars .. + DOUBLE PRECISION BR +* .. +* .. Executable Statements .. +* + IF( R.NE.ZERO ) THEN + BR = B * R + if( BR.NE.ZERO ) THEN + DLADIV2 = (A + BR) * T + ELSE + DLADIV2 = A * T + (B * T) * R + END IF + ELSE + DLADIV2 = (A + D * (B / C)) * T + END IF +* + RETURN +* +* End of DLADIV12 +* + END diff --git a/SRC/sladiv.f b/SRC/sladiv.f index eace5c7..6d26da2 100644 --- a/SRC/sladiv.f +++ b/SRC/sladiv.f @@ -36,8 +36,9 @@ *> p + i*q = --------- *> c + i*d *> -*> The algorithm is due to Robert L. Smith and can be found -*> in D. Knuth, The art of Computer Programming, Vol.2, p.195 +*> The algorithm is due to Michael Baudin and Robert L. Smith +*> and can be found in the paper +*> "A Robust Complex Division in Scilab" *> \endverbatim * * Arguments: @@ -83,17 +84,17 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date September 2012 +*> \date January 2013 * *> \ingroup auxOTHERauxiliary * * ===================================================================== SUBROUTINE SLADIV( A, B, C, D, P, Q ) * -* -- LAPACK auxiliary routine (version 3.4.2) -- +* -- LAPACK auxiliary routine (version 3.5.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* September 2012 +* January 2013 * * .. Scalar Arguments .. REAL A, B, C, D, P, Q @@ -101,24 +102,148 @@ * * ===================================================================== * +* .. Parameters .. + REAL BS + PARAMETER ( BS = 2.0E0 ) + REAL HALF + PARAMETER ( HALF = 0.5E0 ) + REAL TWO + PARAMETER ( TWO = 2.0E0 ) +* * .. Local Scalars .. - REAL E, F + REAL AA, BB, CC, DD, AB, CD, S, OV, UN, BE, EPS +* .. +* .. External Functions .. + REAL SLAMCH + EXTERNAL SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL SLADIV1 * .. * .. Intrinsic Functions .. - INTRINSIC ABS + INTRINSIC ABS, MAX +* .. +* .. Executable Statements .. +* + AA = A + BB = B + CC = C + DD = D + AB = MAX( ABS(A), ABS(B) ) + CD = MAX( ABS(C), ABS(D) ) + S = 1.0E0 + + OV = SLAMCH( 'Overflow threshold' ) + UN = SLAMCH( 'Safe minimum' ) + EPS = SLAMCH( 'Epsilon' ) + BE = BS / (EPS*EPS) + + IF( AB >= HALF*OV ) THEN + AA = HALF * AA + BB = HALF * BB + S = TWO * S + END IF + IF( CD >= HALF*OV ) THEN + CC = HALF * CC + DD = HALF * DD + S = HALF * S + END IF + IF( AB <= UN*BS/EPS ) THEN + AA = AA * BE + BB = BB * BE + S = S / BE + END IF + IF( CD <= UN*BS/EPS ) THEN + CC = CC * BE + DD = DD * BE + S = S * BE + END IF + IF( ABS( D ).LE.ABS( C ) ) THEN + CALL SLADIV1(AA, BB, CC, DD, P, Q) + ELSE + CALL SLADIV1(BB, AA, DD, CC, P, Q) + Q = -Q + END IF + P = P * S + Q = Q * S +* + RETURN +* +* End of SLADIV +* + END + + + + SUBROUTINE SLADIV1( A, B, C, D, P, Q ) +* +* -- LAPACK auxiliary routine (version 3.5.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* January 2013 +* +* .. Scalar Arguments .. + REAL A, B, C, D, P, Q +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE + PARAMETER ( ONE = 1.0E0 ) +* +* .. Local Scalars .. + REAL R, T +* .. +* .. External Functions .. + REAL SLADIV2 + EXTERNAL SLADIV2 +* .. +* .. Executable Statements .. +* + R = D / C + T = ONE / (C + D * R) + P = SLADIV2(A, B, C, D, R, T) + A = -A + Q = SLADIV2(B, A, C, D, R, T) +* + RETURN +* +* End of SLADIV1 +* + END + + REAL FUNCTION SLADIV2( A, B, C, D, R, T ) +* +* -- LAPACK auxiliary routine (version 3.5.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* January 2013 +* +* .. Scalar Arguments .. + REAL A, B, C, D, R, T +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO + PARAMETER ( ZERO = 0.0E0 ) +* +* .. Local Scalars .. + REAL BR * .. * .. Executable Statements .. * - IF( ABS( D ).LT.ABS( C ) ) THEN - E = D / C - F = C + D*E - P = ( A+B*E ) / F - Q = ( B-A*E ) / F + IF( R.NE.ZERO ) THEN + BR = B * R + if( BR.NE.ZERO ) THEN + SLADIV2 = (A + BR) * T + ELSE + SLADIV2 = A * T + (B * T) * R + END IF ELSE - E = C / D - F = D + C*E - P = ( B+A*E ) / F - Q = ( -A+B*E ) / F + SLADIV2 = (A + D * (B / C)) * T END IF * RETURN -- 2.7.4