allow setting off GSL default error handler

10 files changed, 51 insertions, 100 deletions

diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index 1212968..5548285 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+{-# OPTIONS_GHC -fglasgow-exts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Common
@@ -23,13 +23,6 @@ import Data.List(transpose,intersperse)
 import Data.Typeable
 import Data.Maybe(fromJust)
 
-debug x = trace (show x) x
-
-data Vector t = V { dim  :: Int
-                  , fptr :: ForeignPtr t
-                  , ptr  :: Ptr t
-                  } -- deriving Typeable
-
 ----------------------------------------------------------------------
 instance (Storable a, RealFloat a) => Storable (Complex a) where    --
     alignment x = alignment (realPart x)                            --
@@ -40,6 +33,8 @@ instance (Storable a, RealFloat a) => Storable (Complex a) where    --
     poke p (a :+ b) = pokeArray (castPtr p) [a,b]                   --
 ----------------------------------------------------------------------
 
+debug x = trace (show x) x
+
 on :: (a -> a -> b) -> (t -> a) -> t -> t -> b
 on f g = \x y -> f (g x) (g y)
 
@@ -55,13 +50,12 @@ common f = commonval . map f where
     commonval [a] = Just a
     commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
 
-xor :: Bool -> Bool -> Bool
-xor a b = a && not b || b && not a
-
 (//) :: x -> (x -> y) -> y
 infixl 0 //
 (//) = flip ($)
 
+-- our codes should start from 1024
+
 errorCode :: Int -> String
 errorCode 1000 = "bad size"
 errorCode 1001 = "bad function code"
@@ -71,25 +65,6 @@ errorCode 1004 = "singular"
 errorCode 1005 = "didn't converge"
 errorCode n    = "code "++show n
 
-check :: String -> [Vector a] -> IO Int -> IO ()
-check msg ls f = do
-    err <- f
-    when (err/=0) (error (msg++": "++errorCode err))
-    mapM_ (touchForeignPtr . fptr) ls
-    return ()
-
---class (Storable a, Typeable a) => Field a
---instance (Storable a, Typeable a) => Field a
-
---isReal :: (Data.Typeable.Typeable a) => (t -> a) -> t -> Bool
---isReal w x   = typeOf (undefined :: Double) == typeOf (w x)
-
---isComp :: (Data.Typeable.Typeable a) => (t -> a) -> t -> Bool
---isComp w x = typeOf (undefined :: Complex Double) == typeOf (w x)
-
---scast :: forall a . forall b . (Typeable a, Typeable b) => a -> b
---scast = fromJust . cast
-
 {- | conversion of Haskell functions into function pointers that can be used in the C side
 -}
 foreign import ccall "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) 
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index f2646a4..0d9dc70 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -21,6 +21,28 @@ import Complex
 import Control.Monad(when)
 import Data.List(transpose)
 import Debug.Trace(trace)
+import Foreign.C.String(peekCString)
+import Foreign.C.Types
+
+
+data Vector t = V { dim  :: Int
+                  , fptr :: ForeignPtr t
+                  , ptr  :: Ptr t
+                  }
+
+check :: String -> [Vector a] -> IO Int -> IO ()
+check msg ls f = do
+    err <- f
+    when (err/=0) $ if err > 999 -- FIXME, it should be 1024
+                      then (error (msg++": "++errorCode err))
+                      else do
+                        ps <- gsl_strerror err
+                        s <- peekCString ps
+                        error (msg++": "++s)
+    mapM_ (touchForeignPtr . fptr) ls
+    return ()
+
+foreign import ccall "aux.h gsl_strerror" gsl_strerror :: Int -> IO (Ptr CChar)
 
 type Vc t s = Int -> Ptr t -> s
 -- not yet admitted by my haddock version
@@ -30,8 +52,6 @@ type Vc t s = Int -> Ptr t -> s
 vec :: Vector t -> (Vc t s) -> s
 vec v f = f (dim v) (ptr v)
 
---baseOf v = (v `at` 0)
-
 createVector :: Storable a => Int -> IO (Vector a)
 createVector n = do
     when (n <= 0) $ error ("trying to createVector of dim "++show n)
@@ -86,8 +106,6 @@ infixl 9 @>
 (@>) = at
 
 
-
-
 -- | creates a new Vector by joining a list of Vectors
 join :: Storable t => [Vector t] -> Vector t
 join [] = error "joining zero vectors"
@@ -111,7 +129,6 @@ asReal v = V { dim = 2*dim v, fptr =  castForeignPtr (fptr v), ptr = castPtr (pt
 asComplex :: Vector Double -> Vector (Complex Double)
 asComplex v = V { dim = dim v `div` 2, fptr =  castForeignPtr (fptr v), ptr = castPtr (ptr v) }
 
-
 ----------------------------------------------------------------
 
 liftVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
diff --git a/lib/Data/Packed/Internal/aux.h b/lib/Data/Packed/Internal/aux.h
index d055d35..83111e5 100644
--- a/lib/Data/Packed/Internal/aux.h
+++ b/lib/Data/Packed/Internal/aux.h
@@ -24,3 +24,5 @@ int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r));
 
 int diagR(KRVEC(d),RMAT(r));
 int diagC(KCVEC(d),CMAT(r));
+
+const char * gsl_strerror (const int gsl_errno);
diff --git a/lib/GSL.hs b/lib/GSL.hs
index a2137d3..a42ede8 100644
--- a/lib/GSL.hs
+++ b/lib/GSL.hs
@@ -21,20 +21,20 @@ module LinearAlgebra.Algorithms,
 module LAPACK,
 module GSL.Integration,
 module GSL.Differentiation,
---module GSL.Special,
 module GSL.Fourier,
 module GSL.Polynomials,
 module GSL.Minimization,
 module GSL.Matrix,
 module GSL.Compat,
 module Data.Packed.Plot,
-module Complex
+module Complex,
+
+setErrorHandlerOff
 
 ) where
 
 import Data.Packed.Vector hiding (constant)
 import Data.Packed.Matrix hiding ((><), multiply)
---import Data.Packed.Tensor
 import LinearAlgebra.Algorithms hiding (pnorm)
 import LAPACK
 import GSL.Integration
@@ -47,3 +47,4 @@ import GSL.Matrix
 import GSL.Compat
 import Data.Packed.Plot
 import Complex
+import GSL.Special(setErrorHandlerOff)
diff --git a/lib/GSL/Differentiation.hs b/lib/GSL/Differentiation.hs
index 6fd7565..e8e22d2 100644
--- a/lib/GSL/Differentiation.hs
+++ b/lib/GSL/Differentiation.hs
@@ -23,7 +23,7 @@ module GSL.Differentiation (
 ) where
 
 import Foreign
-import Data.Packed.Internal.Common(mkfun,check,(//))
+import Data.Packed.Internal(mkfun,check,(//))
 
 derivGen ::
     Int                   -- ^ type: 0 central, 1 forward, 2 backward
diff --git a/lib/GSL/Integration.hs b/lib/GSL/Integration.hs
index dd76be8..4152c9d 100644
--- a/lib/GSL/Integration.hs
+++ b/lib/GSL/Integration.hs
@@ -21,7 +21,7 @@ module GSL.Integration (
 ) where
 
 import Foreign
-import Data.Packed.Internal.Common(mkfun,check,(//))
+import Data.Packed.Internal(mkfun,check,(//))
 
 
 --------------------------------------------------------------------
diff --git a/lib/GSL/Special.hs b/lib/GSL/Special.hs
index 9eaf733..bdee1b2 100644
--- a/lib/GSL/Special.hs
+++ b/lib/GSL/Special.hs
@@ -41,6 +41,7 @@ module GSL.Special (
 , module GSL.Special.Synchrotron
 , module GSL.Special.Trig
 , module GSL.Special.Zeta
+, setErrorHandlerOff
 )
 where
 
@@ -74,51 +75,6 @@ import GSL.Special.Trig
 import GSL.Special.Zeta
 
 
--------------------- simple functions --------------------------
-
-{- | The error function (/gsl_sf_erf/), defined as 2\/ \\sqrt \\pi * \int\_0\^t \\exp -t\^2 dt.
-
-@> erf 1.5
-0.9661051464753108@
-
--}
-foreign import ccall "gsl-aux.h gsl_sf_erf" erf :: Double -> Double
-
-{- | The Gaussian probability density function (/gsl_sf_erf_Z/), defined as (1\/\\sqrt\{2\\pi\}) \\exp(-x\^2\/2).
-
->> erf_Z 1.5
->0.12951759566589172
-
--}
-foreign import ccall "gsl-aux.h gsl_sf_erf_Z" erf_Z :: Double -> Double
-
--------------------- functions returning sf_result -------------
-
-{- | The regular cylindrical Bessel function of zeroth order, J_0(x). This is
-   the example in the GSL manual, returning the value of the function and
-   the error term: 
-
-@\> bessel_J0_e 5.0
-(-0.1775967713143383,1.9302109579684196e-16)@
-
--}
-bessel_J0_e :: Double -> (Double,Double)
-bessel_J0_e x = createSFR "bessel_J0_e" (gsl_sf_bessel_J0_e x)
-foreign import ccall "gsl-aux.h gsl_sf_bessel_J0_e" gsl_sf_bessel_J0_e :: Double -> Ptr Double -> IO Int
-
--------------------- functions returning sf_result_e10 -------------
-
-{- | (From the GSL manual) \"This function computes the exponential \exp(x) using the @gsl_sf_result_e10@ type to return a result with extended range. This function may be useful if the value of \exp(x) would overflow the numeric range of double\". 
-
-For example:
-
-@\> exp_e10_e 30.0
-(1.0686474581524432,1.4711818964275088e-14,13)@
-
-@\> exp 30.0
-1.0686474581524463e13@
-
--}
-exp_e10_e :: Double -> (Double,Int,Double)
-exp_e10_e x = createSFR_E10 "exp_e10_e" (gsl_sf_exp_e10_e x)
-foreign import ccall "gsl-aux.h gsl_sf_exp_e10_e" gsl_sf_exp_e10_e :: Double -> Ptr Double -> IO Int
+-- | This action removes the GSL default error handler which aborts the program, so that
+-- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort.
+foreign import ccall "gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO ()
diff --git a/lib/GSL/Special/Internal.hs b/lib/GSL/Special/Internal.hs
index c7455d9..ad1aaa2 100644
--- a/lib/GSL/Special/Internal.hs
+++ b/lib/GSL/Special/Internal.hs
@@ -25,7 +25,7 @@ module GSL.Special.Internal (
 where
 
 import Foreign
-import Data.Packed.Internal.Common(check,(//))
+import Data.Packed.Internal(check,(//))
 
 
 data Precision = PrecDouble | PrecSingle | PrecApprox
@@ -39,10 +39,10 @@ type Gsl_mode_t = Int
 
 ----------------------------------------------------------------
 -- | access to a sf_result
-createSFR :: Storable a => t -> (Ptr a -> IO Int) -> (a, a)
+createSFR :: Storable a => String -> (Ptr a -> IO Int) -> (a, a)
 createSFR s f = unsafePerformIO $ do
     p <- mallocArray 2
-    f p // check "createSFR" []
+    f p // check s []
     [val,err] <- peekArray 2 p
     free p
     return (val,err)
@@ -52,12 +52,12 @@ createSFR s f = unsafePerformIO $ do
 -- the sf_result_e10 contains two doubles and the exponent
 
 -- | acces to sf_result_e10
-createSFR_E10 :: (Storable t2, Storable t3, Storable t1) => t -> (Ptr a -> IO Int) -> (t1, t2, t3)
+createSFR_E10 :: (Storable t2, Storable t3, Storable t1) => String -> (Ptr a -> IO Int) -> (t1, t2, t3)
 createSFR_E10 s f = unsafePerformIO $ do
     let sd = sizeOf (0::Double)
     let si = sizeOf (0::Int)
     p <- mallocBytes (2*sd + si)
-    f p // check "createSFR_E10" []
+    f p // check s []
     val <- peekByteOff p 0
     err <- peekByteOff p sd
     expo <- peekByteOff p (2*sd) 
diff --git a/lib/GSL/gsl-aux.c b/lib/GSL/gsl-aux.c
index fc95e6f..eec3e5f 100644
--- a/lib/GSL/gsl-aux.c
+++ b/lib/GSL/gsl-aux.c
@@ -64,6 +64,11 @@
 #define BAD_FILE 1003
 
 
+void no_abort_on_error() {
+    gsl_set_error_handler_off();
+}
+
+
 int toScalarR(int code, KRVEC(x), RVEC(r)) { 
     REQUIRES(rn==1,BAD_SIZE);
     DEBUGMSG("toScalarR");
diff --git a/lib/GSL/gsl-aux.h b/lib/GSL/gsl-aux.h
index 2be07c1..8f8618a 100644
--- a/lib/GSL/gsl-aux.h
+++ b/lib/GSL/gsl-aux.h
@@ -10,6 +10,8 @@
 #define KCVEC(A) int A##n, const gsl_complex*A##p
 #define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p
 
+void no_abort_on_error();
+
 int toScalarR(int code, KRVEC(x), RVEC(r));
 /* norm2, absdif, maximum, posmax, etc. */
 
@@ -58,10 +60,3 @@ int minimizeWithDeriv(double f(int, double*), void df(int, double*, double*),
                       KRVEC(xi), RMAT(sol));
 
 int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr);
-
-double gsl_sf_erf(double);
-double gsl_sf_erf_Z(double);
-double gsl_sf_gamma(double);
-
-int gsl_sf_bessel_J0_e(double, double*); // hmmm...
-int gsl_sf_exp_e10_e(double, double*);   // HMMMMM...