[BACK]Return to dladiv.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dladiv.f between versions 1.12 and 1.13

version 1.12, 2012/12/14 14:22:32 version 1.13, 2014/01/27 09:24:34
Line 36 Line 36
 *>            p + i*q = ---------  *>            p + i*q = ---------
 *>                       c + i*d  *>                       c + i*d
 *>  *>
 *> The algorithm is due to Robert L. Smith and can be found  *> The algorithm is due to Michael Baudin and Robert L. Smith
 *> in D. Knuth, The art of Computer Programming, Vol.2, p.195  *> and can be found in the paper
   *> "A Robust Complex Division in Scilab"
 *> \endverbatim  *> \endverbatim
 *  *
 *  Arguments:  *  Arguments:
Line 83 Line 84
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver 
 *> \author NAG Ltd.   *> \author NAG Ltd. 
 *  *
 *> \date September 2012  *> \date January 2013
 *  *
 *> \ingroup auxOTHERauxiliary  *> \ingroup auxOTHERauxiliary
 *  *
 *  =====================================================================  *  =====================================================================
       SUBROUTINE DLADIV( A, B, C, D, P, Q )        SUBROUTINE DLADIV( 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,    --  *  -- 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..--
 *     September 2012  *     January 2013
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       DOUBLE PRECISION   A, B, C, D, P, Q        DOUBLE PRECISION   A, B, C, D, P, Q
Line 101 Line 102
 *  *
 *  =====================================================================  *  =====================================================================
 *  *
   *     .. Parameters ..
         DOUBLE PRECISION   BS
         PARAMETER          ( BS = 2.0D0 )
         DOUBLE PRECISION   HALF
         PARAMETER          ( HALF = 0.5D0 )
         DOUBLE PRECISION   TWO
         PARAMETER          ( TWO = 2.0D0 )
   *
 *     .. Local Scalars ..  *     .. Local Scalars ..
       DOUBLE PRECISION   E, F        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 Functions ..
       INTRINSIC          ABS        INTRINSIC          ABS, MAX
 *     ..  *     ..
 *     .. Executable Statements ..  *     .. Executable Statements ..
 *  *
       IF( ABS( D ).LT.ABS( C ) ) THEN        AA = A
          E = D / C        BB = B
          F = C + D*E        CC = C
          P = ( A+B*E ) / F        DD = D
          Q = ( B-A*E ) / F        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        ELSE
          E = C / D           CALL DLADIV1(BB, AA, DD, CC, P, Q)
          F = D + C*E           Q = -Q
          P = ( B+A*E ) / F  
          Q = ( -A+B*E ) / F  
       END IF        END IF
         P = P * S
         Q = Q * S
 *  *
       RETURN        RETURN
 *  *
 *     End of DLADIV  *     End of DLADIV
 *  *
       END        END
   
         
   
         SUBROUTINE DLADIV1( 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 ..
         DOUBLE PRECISION   A, B, C, D, P, Q
   *     ..
   *
   *  =====================================================================
   *
   *     .. Parameters ..
         DOUBLE PRECISION   ONE
         PARAMETER          ( ONE = 1.0D0 )
   *
   *     .. Local Scalars ..
         DOUBLE PRECISION   R, T
   *     ..
   *     .. External Functions ..
         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
   
         DOUBLE PRECISION FUNCTION DLADIV2( 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 ..
         DOUBLE PRECISION   A, B, C, D, R, T
   *     ..
   *
   *  =====================================================================
   *
   *     .. Parameters ..
         DOUBLE PRECISION   ZERO
         PARAMETER          ( ZERO = 0.0D0 )
   *
   *     .. Local Scalars ..
         DOUBLE PRECISION   BR
   *     ..
   *     .. Executable Statements ..
   *
         IF( R.NE.ZERO ) THEN
            BR = B * R
            if( BR.NE.ZERO ) THEN
               DLADIV2 = (A + BR) * T
            ELSE
               DLADIV2 = A * T + (B * T) * R
            END IF
         ELSE
            DLADIV2 = (A + D * (B / C)) * T
         END IF
   *
         RETURN
   *
   *     End of DLADIV12
   *
         END
Removed from v.1.12  
changed lines
  Added in v.1.13

CVSweb interface <joel.bertrand@systella.fr>