1 #ifndef STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_ACTION_HANDLER_HPP 2 #define STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_ACTION_HANDLER_HPP 20 static constexpr
int p_max = 8;
21 static constexpr
int m_max = 55;
22 static constexpr
double tol = 1.1e-16;
25 const std::vector<double> theta_m_single_precision{
26 1.3e-1, 1.0e0, 2.2e0, 3.6e0, 4.9e0, 6.3e0,
27 7.7e0, 9.1e0, 1.1e1, 1.2e1, 1.3e1};
28 const std::vector<double> theta_m_double_precision{
29 2.4e-3, 1.4e-1, 6.4e-1, 1.4e0, 2.4e0, 3.5e0,
30 4.7e0, 6.0e0, 7.2e0, 8.5e0, 9.9e0};
32 double l1norm(
const Eigen::MatrixXd& m) {
33 return m.colwise().lpNorm<1>().maxCoeff();
47 inline Eigen::MatrixXd
action(
const Eigen::MatrixXd& mat,
48 const Eigen::MatrixXd& b,
49 const double& t = 1.0) {
50 Eigen::MatrixXd A = mat;
51 double mu = A.trace() / A.rows();
52 for (
int i = 0; i < A.rows(); ++i) {
59 Eigen::MatrixXd res(A.rows(), b.cols());
61 for (
int col = 0;
col < b.cols(); ++
col) {
63 Eigen::VectorXd B = b.col(
col);
64 Eigen::VectorXd F = B;
65 const auto eta =
std::exp(t * mu / s);
66 for (
int i = 1; i < s + 1; ++i) {
67 auto c1 = B.template lpNorm<Eigen::Infinity>();
69 for (
int j = 1; j < m + 1; ++j) {
70 B = t * A * B / (s * j);
71 auto c2 = B.template lpNorm<Eigen::Infinity>();
73 if (c1 + c2 < tol * F.template lpNorm<Eigen::Infinity>()) {
105 const double& t,
int& m,
int& s) {
106 if (l1norm(mat) < tol || t < tol) {
112 Eigen::MatrixXd a = mat * t;
113 const double theta_m = theta_m_double_precision.back();
117 double ap =
std::pow(l1norm(a), 1.0 / p_max);
matrix_exp_action_handler()
Eigen::MatrixXd action(const Eigen::MatrixXd &mat, const Eigen::MatrixXd &b, const double &t=1.0)
fvar< T > exp(const fvar< T > &x)
Eigen::Matrix< T, Eigen::Dynamic, 1 > col(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, size_t j)
Return the specified column of the specified matrix using start-at-1 indexing.
void set_approximation_parameter(const Eigen::MatrixXd &mat, const double &t, int &m, int &s)
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
fvar< T > log2(const fvar< T > &x)
Return the base two logarithm of the specified argument.
fvar< T > ceil(const fvar< T > &x)