1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_INTEGRATE_ODE_CVODES_HPP 2 #define STAN_MATH_REV_MAT_FUNCTOR_INTEGRATE_ODE_CVODES_HPP 15 #include <cvodes/cvodes.h> 16 #include <sunlinsol/sunlinsol_dense.h> 73 template <
typename F,
typename T_initial,
typename T_param,
typename T_t0,
75 std::vector<std::vector<
77 integrate(
const F& f,
const std::vector<T_initial>& y0,
const T_t0& t0,
78 const std::vector<T_ts>& ts,
const std::vector<T_param>& theta,
79 const std::vector<double>& x,
const std::vector<int>& x_int,
80 std::ostream* msgs,
double relative_tolerance,
81 double absolute_tolerance,
82 long int max_num_steps) {
86 const char* fun =
"integrate_ode_cvodes";
89 const std::vector<double> ts_dbl =
value_of(ts);
99 check_less(fun,
"initial time", t0_dbl, ts_dbl[0]);
100 if (relative_tolerance <= 0)
102 relative_tolerance,
"",
", must be greater than 0");
103 if (absolute_tolerance <= 0)
105 absolute_tolerance,
"",
", must be greater than 0");
106 if (max_num_steps <= 0)
108 "",
", must be greater than 0");
110 const size_t N = y0.size();
111 const size_t M = theta.size();
112 const size_t S = (initial_var::value ? N : 0) + (param_var::value ? M : 0);
115 ode_data cvodes_data(f, y0, theta, x, x_int, msgs);
117 void* cvodes_mem = CVodeCreate(Lmm);
118 if (cvodes_mem ==
nullptr)
119 throw std::runtime_error(
"CVodeCreate failed to allocate memory");
121 const size_t coupled_size = cvodes_data.coupled_ode_.size();
123 std::vector<std::vector<
127 f, y0, theta, t0, ts, x, x_int, msgs, y);
131 cvodes_data.nv_state_),
136 CVodeSetUserData(cvodes_mem, reinterpret_cast<void*>(&cvodes_data)),
147 CVodeSetLinearSolver(cvodes_mem, cvodes_data.LS_, cvodes_data.A_),
148 "CVodeSetLinearSolver");
150 CVodeSetJacFn(cvodes_mem, &ode_data::cv_jacobian_states),
156 CVodeSensInit(cvodes_mem, static_cast<int>(S), CV_STAGGERED,
157 &ode_data::cv_rhs_sens, cvodes_data.nv_state_sens_),
161 "CVodeSensEEtolerances");
164 double t_init = t0_dbl;
165 for (
size_t n = 0; n < ts.size(); ++n) {
166 double t_final = ts_dbl[n];
167 if (t_final != t_init)
173 CVodeGetSens(cvodes_mem, &t_init, cvodes_data.nv_state_sens_),
176 observer(cvodes_data.coupled_state_, t_final);
179 }
catch (
const std::exception&
e) {
180 CVodeFree(&cvodes_mem);
184 CVodeFree(&cvodes_mem);
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
void check_nonzero_size(const char *function, const char *name, const T_y &y)
Check if the specified matrix/vector is of non-zero size.
void cvodes_set_options(void *cvodes_mem, double rel_tol, double abs_tol, long int max_num_steps)
void check_ordered(const char *function, const char *name, const std::vector< T_y > &y)
Check if the specified vector is sorted into strictly increasing order.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Defines a public enum named value which is defined to be false as the primitive scalar types cannot b...
Observer for the coupled states.
CVODES ode data holder object which is used during CVODES integration for CVODES callbacks.
boost::math::tools::promote_args< double, typename scalar_type< T >::type, typename return_type< Types_pack... >::type >::type type
void cvodes_check_flag(int flag, const char *func_name)
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw an invalid_argument exception with a consistently formatted message.
double e()
Return the base of the natural logarithm.
std::vector< std::vector< typename stan::return_type< T_initial, T_param, T_t0, T_ts >::type > > integrate(const F &f, const std::vector< T_initial > &y0, const T_t0 &t0, const std::vector< T_ts > &ts, const std::vector< T_param > &theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs, double relative_tolerance, double absolute_tolerance, long int max_num_steps)
Return the solutions for the specified system of ordinary differential equations given the specified ...
void check_less(const char *function, const char *name, const T_y &y, const T_high &high)
Check if y is strictly less than high.
Integrator interface for CVODES' ODE solvers (Adams & BDF methods).