summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2018-03-19 16:22:20 +0000
committerDominic Steinitz <dominic@steinitz.org>2018-03-19 16:22:20 +0000
commit15a2c146d677dea7c4e0d7a39cbd6c8e9edb658f (patch)
tree791db6e42312a7ec8c40a1a5e00b0dccfa6b206c /packages
parent5a51a987e014066d019473a68c1ceca9e30a348f (diff)
More removal of the Fortran way
Diffstat (limited to 'packages')
-rw-r--r--packages/sundials/src/Main.hs14
-rw-r--r--packages/sundials/src/helpers.c35
-rw-r--r--packages/sundials/src/helpers.h9
3 files changed, 0 insertions, 58 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index 89d6668..ab5b153 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -84,19 +84,6 @@ vectorToC vec len ptr = do
84 ptr' <- newForeignPtr_ ptr 84 ptr' <- newForeignPtr_ ptr
85 V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec 85 V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec
86 86
87-- Provided you always call your function 'multiEq' then we can
88-- probably solve any set of ODEs! But of course we don't want to
89-- follow the Fortran way of interacting with sundials.
90
91-- foreign export ccall multiEq :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> Ptr CLong -> Ptr CDouble -> Ptr CInt -> IO ()
92
93multiEq :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> Ptr CLong -> Ptr CDouble -> Ptr CInt -> IO ()
94multiEq tPtr yPtr yDotPtr iParPtr rParPtr ierPtr = do
95 t <- peek tPtr
96 y <- vectorFromC 1 yPtr
97 vectorToC (V.map realToFrac $ stiffish (realToFrac t) (V.map realToFrac y)) 1 yDotPtr
98 poke ierPtr 0
99
100stiffish :: Double -> V.Vector Double -> V.Vector Double 87stiffish :: Double -> V.Vector Double -> V.Vector Double
101stiffish t v = V.fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] 88stiffish t v = V.fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
102 where 89 where
@@ -161,7 +148,6 @@ solveOdeC fun f0 = unsafePerformIO $ do
161 /* Here we use the C types defined in helpers.h which tie up with */ 148 /* Here we use the C types defined in helpers.h which tie up with */
162 /* the Haskell types defined in Types */ 149 /* the Haskell types defined in Types */
163 flag = ARKodeInit(arkode_mem, NULL, $fun:(int (* funIO) (double t, BarType y[], BarType dydt[], void * params)), T0, y); 150 flag = ARKodeInit(arkode_mem, NULL, $fun:(int (* funIO) (double t, BarType y[], BarType dydt[], void * params)), T0, y);
164 /* flag = ARKodeInit(arkode_mem, NULL, FARKfi, T0, y); */
165 if (check_flag(&flag, "ARKodeInit", 1)) return 1; 151 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
166 152
167 /* Set routines */ 153 /* Set routines */
diff --git a/packages/sundials/src/helpers.c b/packages/sundials/src/helpers.c
index 6162b71..f3fe40b 100644
--- a/packages/sundials/src/helpers.c
+++ b/packages/sundials/src/helpers.c
@@ -49,41 +49,6 @@ int check_flag(void *flagvalue, const char *funcname, int opt)
49 return 0; 49 return 0;
50} 50}
51 51
52/* f routine to compute the ODE RHS function f(t,y). */
53int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
54{
55 realtype *rdata = (realtype *) user_data; /* cast user_data to realtype */
56 realtype lamda = rdata[0]; /* set shortcut for stiffness parameter */
57 realtype u = NV_Ith_S(y,0); /* access current solution value */
58
59 /* fill in the RHS function: "NV_Ith_S" accesses the 0th entry of ydot */
60 NV_Ith_S(ydot,0) = lamda*u + 1.0/(1.0+t*t) - lamda*atan(t);
61
62 return 0; /* return with success */
63}
64
65int FARK_IMP_FUN(realtype *T, realtype *Y, realtype *YDOT,
66 long int *IPAR, realtype *RPAR, int *IER) {
67 multiEq(T, Y, YDOT, IPAR, RPAR, IER);
68 return 0;
69}
70
71/* C interface to user-supplied FORTRAN function FARKIFUN; see
72 farkode.h for further details */
73int FARKfi(realtype t, N_Vector y, N_Vector ydot, void *user_data) {
74
75 int ier;
76 realtype *ydata, *dydata;
77 FARKUserData ARK_userdata;
78 ydata = N_VGetArrayPointer(y);
79 dydata = N_VGetArrayPointer(ydot);
80 ARK_userdata = (FARKUserData) user_data;
81
82 FARK_IMP_FUN(&t, ydata, dydata, ARK_userdata->ipar,
83 ARK_userdata->rpar, &ier);
84 return(ier);
85}
86
87/* Jacobian routine to compute J(t,y) = df/dy. */ 52/* Jacobian routine to compute J(t,y) = df/dy. */
88int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J, 53int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
89 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) 54 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 7f4ba02..69a3dfe 100644
--- a/packages/sundials/src/helpers.h
+++ b/packages/sundials/src/helpers.h
@@ -21,14 +21,5 @@ typedef struct _N_VectorContent_Serial BazType;
21*/ 21*/
22int check_flag(void *flagvalue, const char *funcname, int opt); 22int check_flag(void *flagvalue, const char *funcname, int opt);
23 23
24/* f routine to compute the ODE RHS function f(t,y). */
25int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
26
27int FARKfi(realtype t, N_Vector y, N_Vector ydot, void *user_data);
28
29/* Jacobian routine to compute J(t,y) = df/dy. */
30int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
31 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
32
33/* check the computed solution */ 24/* check the computed solution */
34int check_ans(N_Vector y, realtype t, realtype rtol, realtype atol); 25int check_ans(N_Vector y, realtype t, realtype rtol, realtype atol);