25 #ifndef EIGEN_MATRIX_FUNCTION_ATOMIC
26 #define EIGEN_MATRIX_FUNCTION_ATOMIC
38 template <
typename MatrixType>
43 typedef typename MatrixType::Scalar Scalar;
44 typedef typename MatrixType::Index Index;
45 typedef typename NumTraits<Scalar>::Real RealScalar;
46 typedef typename internal::stem_function<Scalar>::type StemFunction;
47 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
58 MatrixType
compute(
const MatrixType& A);
67 bool taylorConverged(Index s,
const MatrixType& F,
const MatrixType& Fincr,
const MatrixType& P);
79 MatrixType m_Ashifted;
85 template <
typename MatrixType>
90 m_avgEival = A.trace() / Scalar(RealScalar(m_Arows));
91 m_Ashifted = A - m_avgEival * MatrixType::Identity(m_Arows, m_Arows);
93 MatrixType F = m_f(m_avgEival, 0) * MatrixType::Identity(m_Arows, m_Arows);
94 MatrixType P = m_Ashifted;
96 for (Index s = 1; s < 1.1 * m_Arows + 10; s++) {
97 Fincr = m_f(m_avgEival, static_cast<int>(s)) * P;
99 P = Scalar(RealScalar(1.0/(s + 1))) * P * m_Ashifted;
100 if (taylorConverged(s, F, Fincr, P)) {
104 eigen_assert(
"Taylor series does not converge" && 0);
109 template <
typename MatrixType>
112 const MatrixType N = MatrixType::Identity(m_Arows, m_Arows) - m_Ashifted;
113 VectorType e = VectorType::Ones(m_Arows);
114 N.template triangularView<Upper>().solveInPlace(e);
115 m_mu = e.cwiseAbs().maxCoeff();
119 template <
typename MatrixType>
120 bool MatrixFunctionAtomic<MatrixType>::taylorConverged(Index s,
const MatrixType& F,
121 const MatrixType& Fincr,
const MatrixType& P)
123 const Index n = F.rows();
124 const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff();
125 const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff();
126 if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) {
127 RealScalar delta = 0;
128 RealScalar rfactorial = 1;
129 for (Index r = 0; r < n; r++) {
131 for (Index i = 0; i < n; i++)
132 mx = (std::max)(mx, std::abs(m_f(m_Ashifted(i, i) + m_avgEival, static_cast<int>(s+r))));
134 rfactorial *= RealScalar(r);
135 delta = (std::max)(delta, mx / rfactorial);
137 const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff();
138 if (m_mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm)
146 #endif // EIGEN_MATRIX_FUNCTION_ATOMIC