Now calling a part of the Haskell ODE

4 files changed, 66 insertions, 24 deletions

diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index 3d5f941..bab5710 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -55,6 +55,13 @@ vectorToC vec len ptr = do
   ptr' <- newForeignPtr_ ptr
   V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec
 
+foreign export ccall singleEq :: Double -> Double -> IO Double
+
+singleEq :: Double -> Double -> IO Double
+singleEq t u = return $ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t
+  where
+    lamda = -100.0
+
 solve :: (CDouble -> V.Vector CDouble -> V.Vector CDouble) ->
          V.Vector Double ->
          CDouble ->
@@ -89,29 +96,6 @@ solve fun f0 lambda = unsafePerformIO $ do
                          realtype abstol = 1.0e-10;
                          realtype lamda  = -100.0;     /* stiffness parameter */
 
-                         /* Beginning of stolen code from the Fortran interface */
-
-                         N_Vector F2C_ARKODE_vec;
-                         F2C_ARKODE_vec = NULL;
-                         F2C_ARKODE_vec = N_VNewEmpty_Serial(NEQ); /* was *N */
-                         if (F2C_ARKODE_vec == NULL) return 1;
-
-                         /* Check for required vector operations */
-                         if(F2C_ARKODE_vec->ops->nvgetarraypointer == NULL) {
-                           fprintf(stderr, "Error: getarraypointer vector operation is not implemented.\n\n");
-                           return 1;
-                         }
-                         if(F2C_ARKODE_vec->ops->nvsetarraypointer == NULL) {
-                           fprintf(stderr, "Error: setarraypointer vector operation is not implemented.\n\n");
-                           return 1;
-                         }
-                         if(F2C_ARKODE_vec->ops->nvcloneempty == NULL) {
-                           fprintf(stderr, "Error: cloneempty vector operation is not implemented.\n\n");
-                           return 1;
-                         }
-
-                         /* End of stolen code from the Fortran interface */
-                         
                          /* Initial diagnostics output */
                          printf("\nAnalytical ODE test problem:\n");
                          printf("    lamda = %"GSYM"\n",    lamda);
@@ -130,7 +114,7 @@ solve fun f0 lambda = unsafePerformIO $ do
                          /*    right-hand side function in y'=f(t,y), the inital time T0, and */
                          /*    the initial dependent variable vector y.  Note: since this */
                          /*    problem is fully implicit, we set f_E to NULL and f_I to f. */
-                         flag = ARKodeInit(arkode_mem, NULL, f, T0, y);
+                         flag = ARKodeInit(arkode_mem, NULL, FARKfi, T0, y);
                          if (check_flag(&flag, "ARKodeInit", 1)) return 1;
 
                          /* Set routines */
diff --git a/packages/sundials/src/helpers.c b/packages/sundials/src/helpers.c
index d51310c..eab5ac9 100644
--- a/packages/sundials/src/helpers.c
+++ b/packages/sundials/src/helpers.c
@@ -8,6 +8,12 @@
 #include <sundials/sundials_types.h>       /* definition of type realtype          */
 #include <sundials/sundials_math.h>
 
+#include "farkode.h"
+
+#include <HsFFI.h>
+#include "Main_stub.h"
+
+
 /* Check function return value...
     opt == 0 means SUNDIALS function allocates memory so check if
              returned NULL pointer
@@ -56,6 +62,31 @@ int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
   return 0;                                   /* return with success */
 }
 
+int FARK_IMP_FUN(realtype *T, realtype *Y, realtype *YDOT,
+                 long int *IPAR, realtype *RPAR, int *IER) {
+  realtype t = *T;
+  realtype u = Y[0];
+  realtype lamda = -100.0;
+  YDOT[0] = singleEq(t, u);
+  return 0;
+}
+
+/* C interface to user-supplied FORTRAN function FARKIFUN; see
+   farkode.h for further details */
+int FARKfi(realtype t, N_Vector y, N_Vector ydot, void *user_data) {
+
+  int ier;
+  realtype *ydata, *dydata;
+  FARKUserData ARK_userdata;
+  ydata  = N_VGetArrayPointer(y);
+  dydata = N_VGetArrayPointer(ydot);
+  ARK_userdata = (FARKUserData) user_data;
+
+  FARK_IMP_FUN(&t, ydata, dydata, ARK_userdata->ipar,
+               ARK_userdata->rpar, &ier);
+  return(ier);
+}
+
 /* Jacobian routine to compute J(t,y) = df/dy. */
 int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
         void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
diff --git a/packages/sundials/src/helpers.h b/packages/sundials/src/helpers.h
index 87b0c81..3b50163 100644
--- a/packages/sundials/src/helpers.h
+++ b/packages/sundials/src/helpers.h
@@ -21,6 +21,8 @@ int check_flag(void *flagvalue, const char *funcname, int opt);
 /* f routine to compute the ODE RHS function f(t,y). */
 int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
 
+int FARKfi(realtype t, N_Vector y, N_Vector ydot, void *user_data);
+
 /* Jacobian routine to compute J(t,y) = df/dy. */
 int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
         void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
diff --git a/packages/sundials/sundials.cabal b/packages/sundials/sundials.cabal
index a68f33a..9b2db12 100644
--- a/packages/sundials/sundials.cabal
+++ b/packages/sundials/sundials.cabal
@@ -18,6 +18,16 @@ cabal-version:       >=1.10
 
 extra-source-files:  src/helpers.c, src/helpers.h
 
+executable ExportAndImport.dylib
+  main-is: ExportAndImport.hs
+  other-extensions:    ForeignFunctionInterface
+  build-depends:       base >=4.10
+  hs-source-dirs:      src
+  default-language:    Haskell2010
+  include-dirs:        src
+  ghc-options:         -O2 -shared -fPIC -dynamic
+  extra-libraries:     HSrts-ghc8.2.2
+
 executable sundials
   main-is:             Main.hs
   -- other-modules:       
@@ -31,3 +41,18 @@ executable sundials
   default-language:    Haskell2010
   extra-libraries:     sundials_arkode
   C-sources:           src/helpers.c src/helpers.h
+
+executable hsundials.dylib
+  main-is:             Main.hs
+  -- other-modules:       
+  other-extensions:    QuasiQuotes, TemplateHaskell, MultiWayIf, OverloadedStrings
+  build-depends:       base >=4.10 && <4.11,
+                       inline-c >=0.6 && <0.7,
+                       vector >=0.12 && <0.13,
+                       template-haskell >=2.12 && <2.13,
+                       containers >=0.5 && <0.6
+  hs-source-dirs:      src
+  ghc-options:         -O2 -shared -fPIC -dynamic
+  default-language:    Haskell2010
+  extra-libraries:     sundials_arkode, HSrts-ghc8.2.2
+  C-sources:           src/helpers.c src/helpers.h