dladiv()

subroutine dladiv	(	double precision	a,
		double precision	b,
		double precision	c,
		double precision	d,
		double precision	p,
		double precision	q
	)

DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.

Download DLADIV + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 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"

Parameters

[in]	A	A is DOUBLE PRECISION
[in]	B	B is DOUBLE PRECISION
[in]	C	C is DOUBLE PRECISION
[in]	D	D is DOUBLE PRECISION The scalars a, b, c, and d in the above expression.
[out]	P	P is DOUBLE PRECISION
[out]	Q	Q is DOUBLE PRECISION The scalars p and q in the above expression.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 90 of file dladiv.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   A, B, C, D, P, Q
*     ..
*
*  =====================================================================
*
*     .. 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   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, 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.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
*

Here is the call graph for this function:

Here is the caller graph for this function: