1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_INTEGRATE_ODE_BDF_HPP 2 #define STAN_MATH_REV_MAT_FUNCTOR_INTEGRATE_ODE_BDF_HPP 14 #include <cvodes/cvodes.h> 15 #include <cvodes/cvodes_band.h> 16 #include <cvodes/cvodes_dense.h> 17 #include <nvector/nvector_serial.h> 35 N_Vector* cvodes_state_sens,
36 void* cvodes_mem,
size_t S) {
37 N_VDestroy_Serial(cvodes_state);
38 if (cvodes_state_sens != NULL)
39 N_VDestroyVectorArray_Serial(cvodes_state_sens, S);
40 CVodeFree(&cvodes_mem);
80 template <
typename F,
typename T_initial,
typename T_param>
84 const std::vector<T_initial>& y0,
86 const std::vector<double>& ts,
87 const std::vector<T_param>& theta,
88 const std::vector<double>& x,
89 const std::vector<int>& x_int,
90 std::ostream* msgs = 0,
91 double relative_tolerance = 1
e-10,
92 double absolute_tolerance = 1
e-10,
93 long int max_num_steps = 1e8) {
100 check_finite(
"integrate_ode_bdf",
"parameter vector", theta);
101 check_finite(
"integrate_ode_bdf",
"continuous data", x);
105 check_less(
"integrate_ode_bdf",
"initial time", t0, ts[0]);
106 if (relative_tolerance <= 0)
108 "relative_tolerance,", relative_tolerance,
109 "",
", must be greater than 0");
110 if (absolute_tolerance <= 0)
112 "absolute_tolerance,", absolute_tolerance,
113 "",
", must be greater than 0");
114 if (max_num_steps <= 0)
116 "max_num_steps,", max_num_steps,
117 "",
", must be greater than 0");
119 const size_t N = y0.size();
120 const size_t M = theta.size();
122 const size_t S = (initial_var::value ? N : 0)
123 + (param_var::value ? M : 0);
124 const size_t size = N * (S + 1);
125 std::vector<double> state(
value_of(y0));
126 N_Vector cvodes_state(N_VMake_Serial(N, &state[0]));
127 N_Vector* cvodes_state_sens = NULL;
130 ode_data cvodes_data(f, y0, theta, x, x_int, msgs);
132 void* cvodes_mem = CVodeCreate(CV_BDF, CV_NEWTON);
133 if (cvodes_mem == NULL)
134 throw std::runtime_error(
"CVodeCreate failed to allocate memory");
136 std::vector<std::vector<double> >
137 y_coupled(ts.size(), std::vector<double>(
size, 0));
146 reinterpret_cast<void*>(&cvodes_data)),
150 relative_tolerance, absolute_tolerance,
157 &ode_data::dense_jacobian),
158 "CVDlsSetDenseJacFn");
162 cvodes_state_sens = N_VCloneVectorArray_Serial(S, cvodes_state);
163 for (
size_t s = 0; s < S; s++)
164 N_VConst(RCONST(0.0), cvodes_state_sens[s]);
168 if (initial_var::value) {
169 for (
size_t n = 0; n < N; n++)
170 NV_Ith_S(cvodes_state_sens[n], n) = 1.0;
174 &ode_data::ode_rhs_sens,
179 "CVodeSensEEtolerances");
183 for (
size_t n = 0; n < ts.size(); ++n) {
184 double t_final = ts[n];
185 if (t_final != t_init)
189 std::copy(state.begin(), state.end(), y_coupled[n].begin());
194 for (
size_t s = 0; s < S; s++)
195 std::copy(NV_DATA_S(cvodes_state_sens[s]),
196 NV_DATA_S(cvodes_state_sens[s]) + N,
197 y_coupled[n].begin() + N + s * N);
201 }
catch (
const std::exception&
e) {
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
std::vector< std::vector< typename stan::return_type< T_initial, T_param >::type > > integrate_ode_bdf(const F &f, const std::vector< T_initial > &y0, double t0, const std::vector< double > &ts, const std::vector< T_param > &theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs=0, double relative_tolerance=1e-10, double absolute_tolerance=1e-10, long int max_num_steps=1e8)
Return the solutions for the specified system of ordinary differential equations given the specified ...
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.
Metaprogram to calculate the base scalar return type resulting from promoting all the scalar types of...
void free_cvodes_memory(N_Vector &cvodes_state, N_Vector *cvodes_state_sens, void *cvodes_mem, size_t S)
Free memory allocated for CVODES state, sensitivity, and general memory.
CVODES ode data holder object which is used during CVODES integration for CVODES callbacks.
void cvodes_check_flag(int flag, const std::string &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.
std::vector< std::vector< typename stan::return_type< T_initial, T_param >::type > > decouple_ode_states(const std::vector< std::vector< double > > &y, const std::vector< T_initial > &y0, const std::vector< T_param > &theta)
Takes sensitivity output from integrators and returns results in precomputed_gradients format...
double e()
Return the base of the natural logarithm.
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
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.