diff options
author | Dominic Steinitz <dominic@steinitz.org> | 2018-03-19 16:22:20 +0000 |
---|---|---|
committer | Dominic Steinitz <dominic@steinitz.org> | 2018-03-19 16:22:20 +0000 |
commit | 15a2c146d677dea7c4e0d7a39cbd6c8e9edb658f (patch) | |
tree | 791db6e42312a7ec8c40a1a5e00b0dccfa6b206c /packages | |
parent | 5a51a987e014066d019473a68c1ceca9e30a348f (diff) |
More removal of the Fortran way
Diffstat (limited to 'packages')
-rw-r--r-- | packages/sundials/src/Main.hs | 14 | ||||
-rw-r--r-- | packages/sundials/src/helpers.c | 35 | ||||
-rw-r--r-- | packages/sundials/src/helpers.h | 9 |
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 | |||
93 | multiEq :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> Ptr CLong -> Ptr CDouble -> Ptr CInt -> IO () | ||
94 | multiEq 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 | |||
100 | stiffish :: Double -> V.Vector Double -> V.Vector Double | 87 | stiffish :: Double -> V.Vector Double -> V.Vector Double |
101 | stiffish t v = V.fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] | 88 | stiffish 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). */ | ||
53 | int 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 | |||
65 | int 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 */ | ||
73 | int 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. */ |
88 | int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J, | 53 | int 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 | */ |
22 | int check_flag(void *flagvalue, const char *funcname, int opt); | 22 | int check_flag(void *flagvalue, const char *funcname, int opt); |
23 | 23 | ||
24 | /* f routine to compute the ODE RHS function f(t,y). */ | ||
25 | int f(realtype t, N_Vector y, N_Vector ydot, void *user_data); | ||
26 | |||
27 | int FARKfi(realtype t, N_Vector y, N_Vector ydot, void *user_data); | ||
28 | |||
29 | /* Jacobian routine to compute J(t,y) = df/dy. */ | ||
30 | int 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 */ |
34 | int check_ans(N_Vector y, realtype t, realtype rtol, realtype atol); | 25 | int check_ans(N_Vector y, realtype t, realtype rtol, realtype atol); |