From 2cb51cd3cbc111a0251560516e772b6cece5446d Mon Sep 17 00:00:00 2001 From: Dominic Steinitz Date: Thu, 15 Mar 2018 16:31:00 +0000 Subject: Now calling a part of the Haskell ODE --- packages/sundials/src/Main.hs | 32 ++++++++------------------------ packages/sundials/src/helpers.c | 31 +++++++++++++++++++++++++++++++ packages/sundials/src/helpers.h | 2 ++ packages/sundials/sundials.cabal | 25 +++++++++++++++++++++++++ 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 /* definition of type realtype */ #include +#include "farkode.h" + +#include +#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 -- cgit v1.2.3