5 files changed, 169 insertions, 14 deletions

diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index 8b713c6..a8d418b 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -1,5 +1,3 @@
-{-# OPTIONS_GHC -Wall #-}
-
 {-# LANGUAGE QuasiQuotes #-}
 {-# LANGUAGE TemplateHaskell #-}
 {-# LANGUAGE MultiWayIf #-}
@@ -141,9 +139,9 @@ import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><),
                                                 subMatrix, rows, cols, toLists,
                                                 size, subVector)
 
-import qualified Types as T
-import           Arkode
-import qualified Arkode as B
+import qualified Numeric.Sundials.CLangToHaskellTypes as T
+import           Numeric.Sundials.Arkode
+import qualified Numeric.Sundials.Arkode as B
 import qualified Numeric.Sundials.ODEOpts as SO
 
 
@@ -160,7 +158,7 @@ C.include "<arkode/arkode_direct.h>"          -- access to ARKDls interface
 C.include "<sundials/sundials_types.h>"       -- definition of type realtype
 C.include "<sundials/sundials_math.h>"
 C.include "../../../helpers.h"
-C.include "Arkode_hsc.h"
+C.include "Numeric/Sundials/Arkode_hsc.h"
 
 
 type Jacobian = Double -> Vector Double -> Matrix Double
@@ -448,7 +446,7 @@ solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO
   diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
   diagMut <- V.thaw diagnostics
   -- We need the types that sundials expects. These are tied together
-  -- in 'Types'. FIXME: The Haskell type is currently empty!
+  -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty!
   let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
       funIO x y f _ptr = do
         -- Convert the pointer we get from C (y) to a vector, and then
@@ -516,7 +514,7 @@ solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO
                          /* problem as fully implicit and set f_E to NULL and f_I to f.         */
 
                          /* Here we use the C types defined in helpers.h which tie up with */
-                         /* the Haskell types defined in Types                             */
+                         /* the Haskell types defined in CLangToHaskellTypes                             */
                          if ($(int method) < MIN_DIRK_NUM) {
                            flag = ARKodeInit(arkode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), NULL, T0, y);
                            if (check_flag(&flag, "ARKodeInit", 1)) return 1;
diff --git a/packages/sundials/src/Numeric/Sundials/Arkode.hsc b/packages/sundials/src/Numeric/Sundials/Arkode.hsc
new file mode 100644
index 0000000..f5e5dc1
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/Arkode.hsc
@@ -0,0 +1,120 @@
+|module Numeric.Sundials.Arkode where
+
+import Foreign
+import Foreign.C.Types
+
+
+#include <stdio.h>
+#include <sundials/sundials_nvector.h>
+#include <sundials/sundials_matrix.h>
+#include <nvector/nvector_serial.h>
+#include <sunmatrix/sunmatrix_dense.h>
+#include <arkode/arkode.h>
+#include <cvode/cvode.h>
+
+
+#def typedef struct _generic_N_Vector SunVector;
+#def typedef struct _N_VectorContent_Serial SunContent;
+
+#def typedef struct _generic_SUNMatrix SunMatrix;
+#def typedef struct _SUNMatrixContent_Dense SunMatrixContent;
+
+getContentMatrixPtr :: Storable a => Ptr b -> IO a
+getContentMatrixPtr ptr = (#peek SunMatrix, content) ptr
+
+getNRows :: Ptr b -> IO CInt
+getNRows ptr = (#peek SunMatrixContent, M) ptr
+putNRows :: CInt -> Ptr b -> IO ()
+putNRows nr ptr = (#poke SunMatrixContent, M) ptr nr
+
+getNCols :: Ptr b -> IO CInt
+getNCols ptr = (#peek SunMatrixContent, N) ptr
+putNCols :: CInt -> Ptr b -> IO ()
+putNCols nc ptr = (#poke SunMatrixContent, N) ptr nc
+
+getMatrixData :: Storable a => Ptr b -> IO a
+getMatrixData ptr = (#peek SunMatrixContent, data) ptr
+
+getContentPtr :: Storable a => Ptr b -> IO a
+getContentPtr ptr = (#peek SunVector, content) ptr
+
+getData :: Storable a => Ptr b -> IO a
+getData ptr = (#peek SunContent, data) ptr
+
+cV_ADAMS :: Int
+cV_ADAMS = #const CV_ADAMS
+cV_BDF :: Int
+cV_BDF = #const CV_BDF
+
+arkSMax :: Int
+arkSMax = #const ARK_S_MAX
+
+mIN_DIRK_NUM, mAX_DIRK_NUM :: Int
+mIN_DIRK_NUM = #const MIN_DIRK_NUM
+mAX_DIRK_NUM = #const MAX_DIRK_NUM
+
+-- FIXME: We could just use inline-c instead
+
+-- Butcher table accessors -- implicit
+sDIRK_2_1_2 :: Int
+sDIRK_2_1_2 = #const SDIRK_2_1_2
+bILLINGTON_3_3_2 :: Int
+bILLINGTON_3_3_2 = #const BILLINGTON_3_3_2
+tRBDF2_3_3_2 :: Int
+tRBDF2_3_3_2 = #const TRBDF2_3_3_2
+kVAERNO_4_2_3 :: Int
+kVAERNO_4_2_3 = #const KVAERNO_4_2_3
+aRK324L2SA_DIRK_4_2_3 :: Int
+aRK324L2SA_DIRK_4_2_3 = #const ARK324L2SA_DIRK_4_2_3
+cASH_5_2_4 :: Int
+cASH_5_2_4 = #const CASH_5_2_4
+cASH_5_3_4 :: Int
+cASH_5_3_4 = #const CASH_5_3_4
+sDIRK_5_3_4 :: Int
+sDIRK_5_3_4 = #const SDIRK_5_3_4
+kVAERNO_5_3_4 :: Int
+kVAERNO_5_3_4 = #const KVAERNO_5_3_4
+aRK436L2SA_DIRK_6_3_4 :: Int
+aRK436L2SA_DIRK_6_3_4 = #const ARK436L2SA_DIRK_6_3_4
+kVAERNO_7_4_5 :: Int
+kVAERNO_7_4_5 = #const KVAERNO_7_4_5
+aRK548L2SA_DIRK_8_4_5 :: Int
+aRK548L2SA_DIRK_8_4_5 = #const ARK548L2SA_DIRK_8_4_5
+
+-- #define DEFAULT_DIRK_2          SDIRK_2_1_2
+-- #define DEFAULT_DIRK_3          ARK324L2SA_DIRK_4_2_3
+-- #define DEFAULT_DIRK_4          SDIRK_5_3_4
+-- #define DEFAULT_DIRK_5          ARK548L2SA_DIRK_8_4_5
+
+-- Butcher table accessors -- explicit
+hEUN_EULER_2_1_2 :: Int
+hEUN_EULER_2_1_2 = #const HEUN_EULER_2_1_2
+bOGACKI_SHAMPINE_4_2_3 :: Int
+bOGACKI_SHAMPINE_4_2_3 = #const BOGACKI_SHAMPINE_4_2_3
+aRK324L2SA_ERK_4_2_3 :: Int
+aRK324L2SA_ERK_4_2_3 = #const ARK324L2SA_ERK_4_2_3
+zONNEVELD_5_3_4 :: Int
+zONNEVELD_5_3_4 = #const ZONNEVELD_5_3_4
+aRK436L2SA_ERK_6_3_4 :: Int
+aRK436L2SA_ERK_6_3_4 = #const ARK436L2SA_ERK_6_3_4
+sAYFY_ABURUB_6_3_4 :: Int
+sAYFY_ABURUB_6_3_4 = #const SAYFY_ABURUB_6_3_4
+cASH_KARP_6_4_5 :: Int
+cASH_KARP_6_4_5 = #const CASH_KARP_6_4_5
+fEHLBERG_6_4_5 :: Int
+fEHLBERG_6_4_5 = #const FEHLBERG_6_4_5
+dORMAND_PRINCE_7_4_5 :: Int
+dORMAND_PRINCE_7_4_5 = #const DORMAND_PRINCE_7_4_5
+aRK548L2SA_ERK_8_4_5 :: Int
+aRK548L2SA_ERK_8_4_5 = #const ARK548L2SA_ERK_8_4_5
+vERNER_8_5_6 :: Int
+vERNER_8_5_6 = #const VERNER_8_5_6
+fEHLBERG_13_7_8 :: Int
+fEHLBERG_13_7_8 = #const FEHLBERG_13_7_8
+
+-- #define DEFAULT_ERK_2           HEUN_EULER_2_1_2
+-- #define DEFAULT_ERK_3           BOGACKI_SHAMPINE_4_2_3
+-- #define DEFAULT_ERK_4           ZONNEVELD_5_3_4
+-- #define DEFAULT_ERK_5           CASH_KARP_6_4_5
+-- #define DEFAULT_ERK_6           VERNER_8_5_6
+-- #define DEFAULT_ERK_8           FEHLBERG_13_7_8
diff --git a/packages/sundials/src/Numeric/Sundials/CLangToHaskellTypes.hs b/packages/sundials/src/Numeric/Sundials/CLangToHaskellTypes.hs
new file mode 100644
index 0000000..0908cbe
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/CLangToHaskellTypes.hs
@@ -0,0 +1,37 @@
+{-# LANGUAGE QuasiQuotes #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE EmptyDataDecls #-}
+
+module Numeric.Sundials.CLangToHaskellTypes where
+
+import           Foreign.C.Types
+
+import qualified Language.Haskell.TH as TH
+import qualified Language.C.Types as CT
+import qualified Data.Map as Map
+import           Language.C.Inline.Context
+
+import qualified Data.Vector.Storable as V
+
+
+data SunVector
+data SunMatrix = SunMatrix { rows :: CInt
+                           , cols :: CInt
+                           , vals :: V.Vector CDouble
+                           }
+
+-- | This is true only if configured/ built as 64 bits
+type SunIndexType = CLong
+
+sunTypesTable :: Map.Map CT.TypeSpecifier TH.TypeQ
+sunTypesTable = Map.fromList
+  [
+    (CT.TypeName "sunindextype", [t| SunIndexType |] )
+  , (CT.TypeName "SunVector", [t| SunVector |] )
+  , (CT.TypeName "SunMatrix", [t| SunMatrix |] )
+  ]
+
+sunCtx :: Context
+sunCtx = mempty {ctxTypesTable = sunTypesTable}
+
diff --git a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
index 0871f9b..1cd072f 100644
--- a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
@@ -90,8 +90,8 @@ import           Numeric.LinearAlgebra.Devel (createVector)
 import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows,
                                                 cols, toLists, size, reshape)
 
-import qualified Types as T
-import           Arkode (cV_ADAMS, cV_BDF)
+import qualified Numeric.Sundials.CLangToHaskellTypes as T
+import           Numeric.Sundials.Arkode (cV_ADAMS, cV_BDF)
 import           Numeric.Sundials.ODEOpts (ODEOpts(..), Jacobian)
 import qualified Numeric.Sundials.ODEOpts as SO
 
@@ -109,7 +109,7 @@ C.include "<cvode/cvode_direct.h>"        -- access to CVDls interface
 C.include "<sundials/sundials_types.h>"   -- definition of type realtype
 C.include "<sundials/sundials_math.h>"
 C.include "../../../helpers.h"
-C.include "Arkode_hsc.h"
+C.include "Numeric/Sundials/Arkode_hsc.h"
 
 
 -- | Stepping functions
@@ -252,7 +252,7 @@ solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts
   diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
   diagMut <- V.thaw diagnostics
   -- We need the types that sundials expects. These are tied together
-  -- in 'Types'. FIXME: The Haskell type is currently empty!
+  -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty!
   let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
       funIO x y f _ptr = do
         -- Convert the pointer we get from C (y) to a vector, and then
diff --git a/packages/sundials/src/Numeric/Sundials/ODEOpts.hs b/packages/sundials/src/Numeric/Sundials/ODEOpts.hs
index 538b474..56dc12c 100644
--- a/packages/sundials/src/Numeric/Sundials/ODEOpts.hs
+++ b/packages/sundials/src/Numeric/Sundials/ODEOpts.hs
@@ -10,8 +10,8 @@ import qualified Data.Vector.Storable.Mutable as VM
 
 import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix)
 
-import qualified Types as T
-import qualified Arkode as B
+import qualified Numeric.Sundials.CLangToHaskellTypes as T
+import qualified Numeric.Sundials.Arkode as B
 
 type Jacobian = Double -> Vector Double -> Matrix Double