1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_CVODES_ODE_DATA_HPP 2 #define STAN_MATH_REV_MAT_FUNCTOR_CVODES_ODE_DATA_HPP 7 #include <cvodes/cvodes.h> 8 #include <sunmatrix/sunmatrix_dense.h> 9 #include <sunlinsol/sunlinsol_dense.h> 10 #include <nvector/nvector_serial.h> 26 template <
typename F,
typename T_initial,
typename T_param>
29 const std::vector<T_initial>& y0_;
30 const std::vector<T_param>& theta_;
31 const std::vector<double> theta_dbl_;
34 const std::vector<double>& x_;
35 const std::vector<int>& x_int_;
74 const std::vector<T_param>& theta,
75 const std::vector<double>& x,
const std::vector<int>& x_int,
86 S_((initial_var::value ? N_ : 0) + (param_var::value ? M_ : 0)),
87 coupled_ode_(f, y0, theta, x, x_int, msgs),
88 coupled_state_(coupled_ode_.initial_state()),
89 nv_state_(N_VMake_Serial(N_, &coupled_state_[0])),
90 nv_state_sens_(nullptr),
91 A_(SUNDenseMatrix(N_, N_)),
92 LS_(SUNDenseLinearSolver(nv_state_, A_)) {
94 nv_state_sens_ = N_VCloneVectorArrayEmpty_Serial(S_, nv_state_);
95 for (std::size_t i = 0; i < S_; i++) {
96 NV_DATA_S(nv_state_sens_[i]) = &coupled_state_[N_] + i * N_;
104 N_VDestroy_Serial(nv_state_);
106 N_VDestroyVectorArray_Serial(nv_state_sens_, S_);
113 static int cv_rhs(realtype t, N_Vector y, N_Vector ydot,
void* user_data) {
114 const ode_data* explicit_ode =
static_cast<const ode_data*
>(user_data);
115 explicit_ode->rhs(t, NV_DATA_S(y), NV_DATA_S(ydot));
123 static int cv_rhs_sens(
int Ns, realtype t, N_Vector y, N_Vector ydot,
124 N_Vector* yS, N_Vector* ySdot,
void* user_data,
125 N_Vector tmp1, N_Vector tmp2) {
126 const ode_data* explicit_ode =
static_cast<const ode_data*
>(user_data);
127 explicit_ode->rhs_sens(t, NV_DATA_S(y), yS, ySdot);
138 SUNMatrix J,
void* user_data, N_Vector tmp1,
139 N_Vector tmp2, N_Vector tmp3) {
140 const ode_data* explicit_ode =
static_cast<const ode_data*
>(user_data);
141 return explicit_ode->jacobian_states(t, NV_DATA_S(y), J);
149 inline void rhs(
double t,
const double y[],
double dy_dt[])
const {
150 const std::vector<double> y_vec(y, y + N_);
151 const std::vector<double>& dy_dt_vec
152 = f_(t, y_vec, theta_dbl_, x_, x_int_, msgs_);
155 std::move(dy_dt_vec.begin(), dy_dt_vec.end(), dy_dt);
165 inline int jacobian_states(
double t,
const double y[], SUNMatrix J)
const {
167 const std::vector<var> y_vec_var(y, y + N_);
170 std::vector<double>&& jacobian_y = std::vector<double>(ode_jacobian.
size());
172 std::move(jacobian_y.begin() + N_, jacobian_y.end(), SM_DATA_D(J));
183 inline void rhs_sens(
double t,
const double y[], N_Vector* yS,
184 N_Vector* ySdot)
const {
185 std::vector<double> z(coupled_state_.size());
186 std::vector<double>&& dz_dt = std::vector<double>(coupled_state_.size());
188 for (std::size_t s = 0; s < S_; s++)
189 std::copy(NV_DATA_S(yS[s]), NV_DATA_S(yS[s]) + N_,
190 z.begin() + (s + 1) * N_);
192 for (std::size_t s = 0; s < S_; s++)
193 std::move(dz_dt.begin() + (s + 1) * N_, dz_dt.begin() + (s + 2) * N_,
194 NV_DATA_S(ySdot[s]));
N_Vector * nv_state_sens_
The coupled_ode_system template specialization for unknown initial values and known parameters...
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...
static int cv_rhs_sens(int Ns, realtype t, N_Vector y, N_Vector ydot, N_Vector *yS, N_Vector *ySdot, void *user_data, N_Vector tmp1, N_Vector tmp2)
Implements the function of type CVSensRhsFn which is the RHS of the sensitivity ODE system...
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
std::vector< double > coupled_state_
static int cv_rhs(realtype t, N_Vector y, N_Vector ydot, void *user_data)
Implements the function of type CVRhsFn which is the user-defined ODE RHS passed to CVODES...
CVODES ode data holder object which is used during CVODES integration for CVODES callbacks.
const coupled_ode_system< F, T_initial, T_param > coupled_ode_
cvodes_ode_data(const F &f, const std::vector< T_initial > &y0, const std::vector< T_param > &theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs)
Construct CVODES ode data object to enable callbacks from CVODES during ODE integration.
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
static void recover_memory_nested()
Recover only the memory used for the top nested call.
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
static int cv_jacobian_states(realtype t, N_Vector y, N_Vector fy, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Implements the function of type CVDlsJacFn which is the user-defined callback for CVODES to calculate...
static void start_nested()
Record the current position so that recover_memory_nested() can find it.
size_t size() const
Returns the size of the coupled system.
const kernel_cl< in_buffer, out_buffer, int, int > copy("copy", {indexing_helpers, copy_kernel_code})
See the docs for copy() .