sladiv()

subroutine sladiv	(	real	a,
		real	b,
		real	c,
		real	d,
		real	p,
		real	q
	)

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

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

Purpose:

 SLADIV 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 REAL
[in]	B	B is REAL
[in]	C	C is REAL
[in]	D	D is REAL The scalars a, b, c, and d in the above expression.
[out]	P	P is REAL
[out]	Q	Q is REAL 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 sladiv.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 ..
      REAL               A, B, C, D, P, Q
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               BS
      parameter( bs = 2.0e0 )
      REAL               HALF
      parameter( half = 0.5e0 )
      REAL               TWO
      parameter( two = 2.0e0 )
*
*     .. Local Scalars ..
      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, 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
*

Here is the call graph for this function:

Here is the caller graph for this function: