Commit Victor Liu's suggestion about making the complex division routine more robust

tests and builds seems fine.

Message sent on Oct 18th
I just saw a paper on ArXiv about making the complex division routine more robust:
http://arxiv.org/abs/1210.4539
Second author is actually the original inventor of the current algorithm in Lapack.
I have attached my modified DLADIV routine, which passes all the tests in the build process.

1 files changed, 141 insertions, 16 deletions

diff --git a/SRC/sladiv.f b/SRC/sladiv.f
index eace5c72..6d26da20 100644
--- a/SRC/sladiv.f
+++ b/SRC/sladiv.f
@@ -36,8 +36,9 @@
 *>            p + i*q = ---------
 *>                       c + i*d
 *>
-*> The algorithm is due to Robert L. Smith and can be found
-*> in D. Knuth, The art of Computer Programming, Vol.2, p.195
+*> 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:
@@ -83,17 +84,17 @@
 *> \author Univ. of Colorado Denver 
 *> \author NAG Ltd. 
 *
-*> \date September 2012
+*> \date January 2013
 *
 *> \ingroup auxOTHERauxiliary
 *
 *  =====================================================================
       SUBROUTINE SLADIV( 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,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     September 2012
+*     January 2013
 *
 *     .. Scalar Arguments ..
       REAL               A, B, C, D, P, Q
@@ -101,24 +102,148 @@
 *
 *  =====================================================================
 *
+*     .. Parameters ..
+      REAL               BS
+      PARAMETER          ( BS = 2.0E0 )
+      REAL               HALF
+      PARAMETER          ( HALF = 0.5E0 )
+      REAL               TWO
+      PARAMETER          ( TWO = 2.0E0 )
+*
 *     .. Local Scalars ..
-      REAL               E, F
+      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
+      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( ABS( D ).LT.ABS( C ) ) THEN
-         E = D / C
-         F = C + D*E
-         P = ( A+B*E ) / F
-         Q = ( B-A*E ) / F
+      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
-         E = C / D
-         F = D + C*E
-         P = ( B+A*E ) / F
-         Q = ( -A+B*E ) / F
+         SLADIV2 = (A + D * (B / C)) * T
       END IF
 *
       RETURN