version 1.1.1.1, 2010/01/26 15:22:45
|
version 1.17, 2017/06/17 10:53:52
|
Line 1
|
Line 1
|
|
*> \brief \b DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DLADIV + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dladiv.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dladiv.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dladiv.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DLADIV( A, B, C, D, P, Q ) |
|
* |
|
* .. Scalar Arguments .. |
|
* DOUBLE PRECISION A, B, C, D, P, Q |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DLADIV performs complex division in real arithmetic |
|
*> |
|
*> a + i*b |
|
*> p + i*q = --------- |
|
*> c + i*d |
|
*> |
|
*> 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: |
|
* ========== |
|
* |
|
*> \param[in] A |
|
*> \verbatim |
|
*> A is DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[in] B |
|
*> \verbatim |
|
*> B is DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[in] C |
|
*> \verbatim |
|
*> C is DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[in] D |
|
*> \verbatim |
|
*> D is DOUBLE PRECISION |
|
*> The scalars a, b, c, and d in the above expression. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] P |
|
*> \verbatim |
|
*> P is DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[out] Q |
|
*> \verbatim |
|
*> Q is DOUBLE PRECISION |
|
*> The scalars p and q in the above expression. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date January 2013 |
|
* |
|
*> \ingroup doubleOTHERauxiliary |
|
* |
|
* ===================================================================== |
SUBROUTINE DLADIV( A, B, C, D, P, Q ) |
SUBROUTINE DLADIV( A, B, C, D, P, Q ) |
* |
* |
* -- LAPACK auxiliary routine (version 3.2) -- |
* -- LAPACK auxiliary routine (version 3.7.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..-- |
* November 2006 |
* January 2013 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
DOUBLE PRECISION A, B, C, D, P, Q |
DOUBLE PRECISION A, B, C, D, P, Q |
* .. |
* .. |
* |
* |
* Purpose |
* ===================================================================== |
* ======= |
|
* |
* |
* DLADIV performs complex division in real arithmetic |
* .. Parameters .. |
|
DOUBLE PRECISION BS |
|
PARAMETER ( BS = 2.0D0 ) |
|
DOUBLE PRECISION HALF |
|
PARAMETER ( HALF = 0.5D0 ) |
|
DOUBLE PRECISION TWO |
|
PARAMETER ( TWO = 2.0D0 ) |
* |
* |
* a + i*b |
* .. Local Scalars .. |
* p + i*q = --------- |
DOUBLE PRECISION AA, BB, CC, DD, AB, CD, S, OV, UN, BE, EPS |
* c + i*d |
* .. |
* |
* .. External Functions .. |
* The algorithm is due to Robert L. Smith and can be found |
DOUBLE PRECISION DLAMCH |
* in D. Knuth, The art of Computer Programming, Vol.2, p.195 |
EXTERNAL DLAMCH |
* |
* .. |
* Arguments |
* .. External Subroutines .. |
* ========= |
EXTERNAL DLADIV1 |
* |
* .. |
* A (input) DOUBLE PRECISION |
* .. Intrinsic Functions .. |
* B (input) DOUBLE PRECISION |
INTRINSIC ABS, MAX |
* C (input) DOUBLE PRECISION |
* .. |
* D (input) DOUBLE PRECISION |
* .. Executable Statements .. |
* The scalars a, b, c, and d in the above expression. |
* |
* |
AA = A |
* P (output) DOUBLE PRECISION |
BB = B |
* Q (output) DOUBLE PRECISION |
CC = C |
* The scalars p and q in the above expression. |
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 |
|
CALL DLADIV1(BB, AA, DD, CC, P, Q) |
|
Q = -Q |
|
END IF |
|
P = P * S |
|
Q = Q * S |
|
* |
|
RETURN |
|
* |
|
* End of DLADIV |
|
* |
|
END |
|
|
|
*> \ingroup doubleOTHERauxiliary |
|
|
|
|
|
SUBROUTINE DLADIV1( A, B, C, D, P, Q ) |
|
* |
|
* -- LAPACK auxiliary routine (version 3.7.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 .. |
* .. Local Scalars .. |
DOUBLE PRECISION E, F |
DOUBLE PRECISION R, T |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. External Functions .. |
INTRINSIC ABS |
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 |
|
|
|
*> \ingroup doubleOTHERauxiliary |
|
|
|
DOUBLE PRECISION FUNCTION DLADIV2( A, B, C, D, R, T ) |
|
* |
|
* -- LAPACK auxiliary routine (version 3.7.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 .. |
* .. Executable Statements .. |
* |
* |
IF( ABS( D ).LT.ABS( C ) ) THEN |
IF( R.NE.ZERO ) THEN |
E = D / C |
BR = B * R |
F = C + D*E |
IF( BR.NE.ZERO ) THEN |
P = ( A+B*E ) / F |
DLADIV2 = (A + BR) * T |
Q = ( B-A*E ) / F |
ELSE |
|
DLADIV2 = A * T + (B * T) * R |
|
END IF |
ELSE |
ELSE |
E = C / D |
DLADIV2 = (A + D * (B / C)) * T |
F = D + C*E |
|
P = ( B+A*E ) / F |
|
Q = ( -A+B*E ) / F |
|
END IF |
END IF |
* |
* |
RETURN |
RETURN |
* |
* |
* End of DLADIV |
* End of DLADIV12 |
* |
* |
END |
END |