From 3adf82b7ee54b5866013449bda49cb545c13ee47 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 5 Feb 2014 20:11:43 +0100 Subject: zero absolute tolerance in the integration routines provisionally changed to 1E-12 --- lib/Numeric/GSL/Integration.hs | 51 ++++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 24 deletions(-) diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs index d1f552d..7c1cdb9 100644 --- a/lib/Numeric/GSL/Integration.hs +++ b/lib/Numeric/GSL/Integration.hs @@ -31,6 +31,8 @@ import Foreign.Storable(peek) import Data.Packed.Internal(check,(//)) import System.IO.Unsafe(unsafePerformIO) +eps = 1e-12 + {- | conversion of Haskell functions into function pointers that can be used in the C side -} foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) @@ -55,7 +57,7 @@ integrateQAGS prec n f a b = unsafePerformIO $ do r <- malloc e <- malloc fp <- mkfun (\x _ -> f x) - c_integrate_qags fp a b prec (fromIntegral n) r e // check "integrate_qags" + c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags" vr <- peek r ve <- peek e let result = (vr,ve) @@ -64,9 +66,9 @@ integrateQAGS prec n f a b = unsafePerformIO $ do freeHaskellFunPtr fp return result -foreign import ccall safe "gsl-aux.h integrate_qags" - c_integrate_qags :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> CInt - -> Ptr Double -> Ptr Double -> IO CInt +foreign import ccall safe "integrate_qags" c_integrate_qags + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt ----------------------------------------------------------------- {- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example: @@ -86,7 +88,7 @@ integrateQNG prec f a b = unsafePerformIO $ do r <- malloc e <- malloc fp <- mkfun (\x _ -> f x) - c_integrate_qng fp a b prec r e // check "integrate_qng" + c_integrate_qng fp a b eps prec r e // check "integrate_qng" vr <- peek r ve <- peek e let result = (vr,ve) @@ -95,9 +97,9 @@ integrateQNG prec f a b = unsafePerformIO $ do freeHaskellFunPtr fp return result -foreign import ccall safe "gsl-aux.h integrate_qng" - c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double - -> Ptr Double -> Ptr Double -> IO CInt +foreign import ccall safe "integrate_qng" c_integrate_qng + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt -------------------------------------------------------------------- {- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS). @@ -118,7 +120,7 @@ integrateQAGI prec n f = unsafePerformIO $ do r <- malloc e <- malloc fp <- mkfun (\x _ -> f x) - c_integrate_qagi fp prec (fromIntegral n) r e // check "integrate_qagi" + c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi" vr <- peek r ve <- peek e let result = (vr,ve) @@ -127,9 +129,9 @@ integrateQAGI prec n f = unsafePerformIO $ do freeHaskellFunPtr fp return result -foreign import ccall safe "gsl-aux.h integrate_qagi" - c_integrate_qagi :: FunPtr (Double-> Ptr() -> Double) -> Double -> CInt - -> Ptr Double -> Ptr Double -> IO CInt +foreign import ccall safe "integrate_qagi" c_integrate_qagi + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> CInt -> Ptr Double -> Ptr Double -> IO CInt -------------------------------------------------------------------- {- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf). @@ -151,7 +153,7 @@ integrateQAGIU prec n f a = unsafePerformIO $ do r <- malloc e <- malloc fp <- mkfun (\x _ -> f x) - c_integrate_qagiu fp a prec (fromIntegral n) r e // check "integrate_qagiu" + c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu" vr <- peek r ve <- peek e let result = (vr,ve) @@ -160,9 +162,9 @@ integrateQAGIU prec n f a = unsafePerformIO $ do freeHaskellFunPtr fp return result -foreign import ccall safe "gsl-aux.h integrate_qagiu" - c_integrate_qagiu :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt - -> Ptr Double -> Ptr Double -> IO CInt +foreign import ccall safe "integrate_qagiu" c_integrate_qagiu + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt -------------------------------------------------------------------- {- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b). @@ -184,7 +186,7 @@ integrateQAGIL prec n f b = unsafePerformIO $ do r <- malloc e <- malloc fp <- mkfun (\x _ -> f x) - c_integrate_qagil fp b prec (fromIntegral n) r e // check "integrate_qagil" + c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil" vr <- peek r ve <- peek e let result = (vr,ve) @@ -193,9 +195,9 @@ integrateQAGIL prec n f b = unsafePerformIO $ do freeHaskellFunPtr fp return result -foreign import ccall safe "gsl-aux.h integrate_qagil" - c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt - -> Ptr Double -> Ptr Double -> IO CInt +foreign import ccall safe "gsl-aux.h integrate_qagil" c_integrate_qagil + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt -------------------------------------------------------------------- @@ -231,7 +233,7 @@ integrateCQUAD prec n f a b = unsafePerformIO $ do e <- malloc neval <- malloc fp <- mkfun (\x _ -> f x) - c_integrate_cquad fp a b prec (fromIntegral n) r e neval // check "integrate_cquad" + c_integrate_cquad fp a b eps prec (fromIntegral n) r e neval // check "integrate_cquad" vr <- peek r ve <- peek e vneval <- peek neval @@ -242,6 +244,7 @@ integrateCQUAD prec n f a b = unsafePerformIO $ do freeHaskellFunPtr fp return result -foreign import ccall safe "gsl-aux.h integrate_cquad" - c_integrate_cquad :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> CInt - -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt +foreign import ccall safe "integrate_cquad" c_integrate_cquad + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt + -- cgit v1.2.3