1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_IDAS_INTEGRATOR_HPP 2 #define STAN_MATH_REV_MAT_FUNCTOR_IDAS_INTEGRATOR_HPP 10 #include <idas/idas.h> 11 #include <sunmatrix/sunmatrix_dense.h> 12 #include <sunlinsol/sunlinsol_dense.h> 13 #include <nvector/nvector_serial.h> 29 const int64_t max_num_steps_;
33 template <
typename Dae>
34 void init_sensitivity(Dae& dae);
64 const double& t0,
const std::vector<double>& ts,
65 std::vector<std::vector<double> >& res_yy);
67 template <
typename Dae>
68 void solve(Dae& dae,
const double& t0,
const std::vector<double>& ts,
69 typename Dae::return_type& res_yy);
82 const int64_t max_num_steps = IDAS_MAX_STEPS)
83 : rtol_(rtol), atol_(atol), max_num_steps_(max_num_steps) {
86 ", must be greater than 0");
89 ", must be less than 1.0E-3");
92 ", must be greater than 0");
93 if (max_num_steps_ <= 0)
95 ", must be greater than 0");
113 template <
typename Dae>
114 typename Dae::return_type
integrate(Dae& dae,
double t0,
115 const std::vector<double>& ts) {
116 using Eigen::Dynamic;
118 using Eigen::MatrixXd;
119 using Eigen::VectorXd;
121 static const char* caller =
"idas_integrator";
126 check_less(caller,
"initial time", t0, ts.front());
128 auto mem = dae.mem();
129 auto yy = dae.nv_yy();
130 auto yp = dae.nv_yp();
131 const size_t n = dae.n();
133 typename Dae::return_type res_yy(
134 ts.size(), std::vector<typename Dae::scalar_type>(n, 0));
136 auto A = SUNDenseMatrix(n, n);
137 auto LS = SUNDenseLinearSolver(yy, A);
147 init_sensitivity(dae);
149 solve(dae, t0, ts, res_yy);
150 }
catch (
const std::exception&
e) {
171 template <
typename Dae>
172 void idas_integrator::init_sensitivity(Dae& dae) {
173 if (Dae::need_sens) {
174 auto mem = dae.mem();
175 auto yys = dae.nv_yys();
176 auto yps = dae.nv_yps();
179 if (Dae::is_var_yy0) {
180 for (
size_t i = 0; i < n; ++i)
181 NV_Ith_S(yys[i], i) = 1.0;
183 if (Dae::is_var_yp0) {
184 for (
size_t i = 0; i < n; ++i)
185 NV_Ith_S(yps[i + n], i) = 1.0;
188 dae.sensitivity_residual(), yys, yps));
203 template <
typename F>
205 const double& t0,
const std::vector<double>& ts,
206 std::vector<std::vector<double> >& res_yy) {
209 auto mem = dae.
mem();
210 auto yy = dae.
nv_yy();
211 auto yp = dae.
nv_yp();
213 std::for_each(ts.begin(), ts.end(), [&](
double t2) {
230 template <
typename Dae>
231 void idas_integrator::solve(Dae& dae,
const double& t0,
232 const std::vector<double>& ts,
233 typename Dae::return_type& res_yy) {
236 auto mem = dae.mem();
237 auto yy = dae.nv_yy();
238 auto yp = dae.nv_yp();
239 auto yys = dae.nv_yys();
240 const auto n = dae.n();
241 const auto ns = dae.ns();
242 auto vars = dae.vars();
244 std::vector<stan::math::var> sol_t(n);
245 std::vector<double> sol_grad(ns);
247 std::for_each(ts.begin(), ts.end(), [&](
double t2) {
250 for (
size_t k = 0; k < n; ++k) {
251 for (
size_t j = 0; j < ns; ++j) {
252 sol_grad[j] = NV_Ith_S(yys[j], k);
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.
Dae::return_type integrate(Dae &dae, double t0, const std::vector< double > &ts)
Return the solutions for the specified DAE given the specified initial state, initial times...
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.
idas_integrator(const double rtol, const double atol, const int64_t max_num_steps=IDAS_MAX_STEPS)
constructor
IDAS DAE system with forward sensitivity calculation.
static constexpr int IDAS_MAX_STEPS
var precomputed_gradients(double value, const std::vector< var > &operands, const std::vector< double > &gradients)
This function returns a var for an expression that has the specified value, vector of operands...
void * mem()
return IDAS memory handle
N_Vector & nv_yp()
return reference to current N_Vector of derivative variable
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.
#define CHECK_IDAS_CALL(call)
double e()
Return the base of the natural logarithm.
const double E
The base of the natural logarithm, .
N_Vector & nv_yy()
return reference to current N_Vector of unknown variable
const std::vector< double > & yy_val()
return reference to current solution vector value
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.
const kernel_cl< in_buffer, out_buffer, int, int > copy("copy", {indexing_helpers, copy_kernel_code})
See the docs for copy() .