From f2b1eae3d71c546abc71e099b4bd86010627f0fb Mon Sep 17 00:00:00 2001
From: Dominic Steinitz <dominic@steinitz.org>
Date: Mon, 12 Mar 2018 09:50:59 +0000
Subject: Now builds with stack and cabal

---
 packages/sundials/src/helpers.c | 36 +++++++++++++++++++++++++++++++++++-
 1 file changed, 35 insertions(+), 1 deletion(-)

(limited to 'packages/sundials/src/helpers.c')

diff --git a/packages/sundials/src/helpers.c b/packages/sundials/src/helpers.c
index 3498119..3e1a9a8 100644
--- a/packages/sundials/src/helpers.c
+++ b/packages/sundials/src/helpers.c
@@ -1,5 +1,12 @@
 #include <stdio.h>
-#include "helpers.h"
+#include <math.h>
+#include <arkode/arkode.h>                 /* prototypes for ARKODE fcts., consts. */
+#include <nvector/nvector_serial.h>        /* serial N_Vector types, fcts., macros */
+#include <sunmatrix/sunmatrix_dense.h>     /* access to dense SUNMatrix            */
+#include <sunlinsol/sunlinsol_dense.h>     /* access to dense SUNLinearSolver      */
+#include <arkode/arkode_direct.h>          /* access to ARKDls interface           */
+#include <sundials/sundials_types.h>       /* definition of type realtype          */
+#include <sundials/sundials_math.h>
 
 /* Check function return value...
     opt == 0 means SUNDIALS function allocates memory so check if
@@ -35,3 +42,30 @@ int check_flag(void *flagvalue, const char *funcname, int opt)
 
   return 0;
 }
+
+/* f routine to compute the ODE RHS function f(t,y). */
+int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
+{
+  realtype *rdata = (realtype *) user_data;   /* cast user_data to realtype */
+  realtype lamda = rdata[0];                  /* set shortcut for stiffness parameter */
+  realtype u = NV_Ith_S(y,0);                 /* access current solution value */
+
+  /* fill in the RHS function: "NV_Ith_S" accesses the 0th entry of ydot */
+  NV_Ith_S(ydot,0) = lamda*u + 1.0/(1.0+t*t) - lamda*atan(t);
+
+  return 0;                                   /* return with success */
+}
+
+/* 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)
+{
+  realtype *rdata = (realtype *) user_data;   /* cast user_data to realtype */
+  realtype lamda = rdata[0];                  /* set shortcut for stiffness parameter */
+  realtype *Jdata = SUNDenseMatrix_Data(J);
+  
+  /* Fill in Jacobian of f: set the first entry of the data array to set the (0,0) entry */
+  Jdata[0] = lamda;
+
+  return 0;                                   /* return with success */
+}
-- 
cgit v1.2.3