version 1.12, 2012/12/14 14:22:32
|
version 1.13, 2014/01/27 09:24:34
|
Line 36
|
Line 36
|
*> p + i*q = --------- |
*> p + i*q = --------- |
*> c + i*d |
*> c + i*d |
*> |
*> |
*> The algorithm is due to Robert L. Smith and can be found |
*> The algorithm is due to Michael Baudin and Robert L. Smith |
*> in D. Knuth, The art of Computer Programming, Vol.2, p.195 |
*> and can be found in the paper |
|
*> "A Robust Complex Division in Scilab" |
*> \endverbatim |
*> \endverbatim |
* |
* |
* Arguments: |
* Arguments: |
Line 83
|
Line 84
|
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date September 2012 |
*> \date January 2013 |
* |
* |
*> \ingroup auxOTHERauxiliary |
*> \ingroup auxOTHERauxiliary |
* |
* |
* ===================================================================== |
* ===================================================================== |
SUBROUTINE DLADIV( A, B, C, D, P, Q ) |
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, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* September 2012 |
* January 2013 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
DOUBLE PRECISION A, B, C, D, P, Q |
DOUBLE PRECISION A, B, C, D, P, Q |
Line 101
|
Line 102
|
* |
* |
* ===================================================================== |
* ===================================================================== |
* |
* |
|
* .. Parameters .. |
|
DOUBLE PRECISION BS |
|
PARAMETER ( BS = 2.0D0 ) |
|
DOUBLE PRECISION HALF |
|
PARAMETER ( HALF = 0.5D0 ) |
|
DOUBLE PRECISION TWO |
|
PARAMETER ( TWO = 2.0D0 ) |
|
* |
* .. Local Scalars .. |
* .. 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 Functions .. |
INTRINSIC ABS |
INTRINSIC ABS, MAX |
* .. |
* .. |
* .. Executable Statements .. |
* .. Executable Statements .. |
* |
* |
IF( ABS( D ).LT.ABS( C ) ) THEN |
AA = A |
E = D / C |
BB = B |
F = C + D*E |
CC = C |
P = ( A+B*E ) / F |
DD = D |
Q = ( B-A*E ) / F |
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 |
ELSE |
E = C / D |
CALL DLADIV1(BB, AA, DD, CC, P, Q) |
F = D + C*E |
Q = -Q |
P = ( B+A*E ) / F |
|
Q = ( -A+B*E ) / F |
|
END IF |
END IF |
|
P = P * S |
|
Q = Q * S |
* |
* |
RETURN |
RETURN |
* |
* |
* End of DLADIV |
* End of DLADIV |
* |
* |
END |
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 |