sladiv.f

Go to the documentation of this file.

*> \brief \b SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download SLADIV + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sladiv.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sladiv.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sladiv.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE SLADIV( A, B, C, D, P, Q )
* 
*       .. Scalar Arguments ..
*       REAL               A, B, C, D, P, Q
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> 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"
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] A
*> \verbatim
*>          A is REAL
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is REAL
*> \endverbatim
*>
*> \param[in] C
*> \verbatim
*>          C is REAL
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is REAL
*>          The scalars a, b, c, and d in the above expression.
*> \endverbatim
*>
*> \param[out] P
*> \verbatim
*>          P is REAL
*> \endverbatim
*>
*> \param[out] Q
*> \verbatim
*>          Q is REAL
*>          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 auxOTHERauxiliary
*
*  =====================================================================
      SUBROUTINE sladiv( 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               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
*
      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( 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
         sladiv2 = (a + d * (b / c)) * t
      END IF
*
      RETURN
*
*     End of SLADIV
*
      END