diff options
Diffstat (limited to 'unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h')
| -rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h | 44 |
1 files changed, 28 insertions, 16 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 6825a7882..f8bbd59c1 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -157,6 +157,24 @@ class MatrixExponential { /** \brief L1 norm of m_M. */ RealScalar m_l1norm; + + /** \brief Scaling operator. */ + struct ScalingOp; +}; + +template <typename MatrixType> +struct MatrixExponential<MatrixType>::ScalingOp +{ + ScalingOp(int squarings) : m_squarings(squarings) { } + + inline const RealScalar operator() (const RealScalar& x) const + { return std::ldexp(x, -m_squarings); } + + inline const ComplexScalar operator() (const ComplexScalar& x) const + { return ComplexScalar(std::ldexp(x.real(), -m_squarings), std::ldexp(x.imag(), -m_squarings)); } + + private: + int m_squarings; }; template <typename MatrixType> @@ -283,17 +301,15 @@ EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade17(const MatrixType template <typename MatrixType> void MatrixExponential<MatrixType>::computeUV(float) { - using std::frexp; - using std::pow; if (m_l1norm < 4.258730016922831e-001) { pade3(m_M); } else if (m_l1norm < 1.880152677804762e+000) { pade5(m_M); } else { const float maxnorm = 3.925724783138660f; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp<ScalingOp, MatrixType>(m_M, m_squarings); pade7(A); } } @@ -301,8 +317,6 @@ void MatrixExponential<MatrixType>::computeUV(float) template <typename MatrixType> void MatrixExponential<MatrixType>::computeUV(double) { - using std::frexp; - using std::pow; if (m_l1norm < 1.495585217958292e-002) { pade3(m_M); } else if (m_l1norm < 2.539398330063230e-001) { @@ -313,9 +327,9 @@ void MatrixExponential<MatrixType>::computeUV(double) pade9(m_M); } else { const double maxnorm = 5.371920351148152; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp<ScalingOp, MatrixType>(m_M, m_squarings); pade13(A); } } @@ -323,8 +337,6 @@ void MatrixExponential<MatrixType>::computeUV(double) template <typename MatrixType> void MatrixExponential<MatrixType>::computeUV(long double) { - using std::frexp; - using std::pow; #if LDBL_MANT_DIG == 53 // double precision computeUV(double()); #elif LDBL_MANT_DIG <= 64 // extended precision @@ -338,9 +350,9 @@ void MatrixExponential<MatrixType>::computeUV(long double) pade9(m_M); } else { const long double maxnorm = 4.0246098906697353063L; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp<ScalingOp, MatrixType>(m_M, m_squarings); pade13(A); } #elif LDBL_MANT_DIG <= 106 // double-double @@ -356,9 +368,9 @@ void MatrixExponential<MatrixType>::computeUV(long double) pade13(m_M); } else { const long double maxnorm = 3.2579440895405400856599663723517L; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp<ScalingOp, MatrixType>(m_M, m_squarings); pade17(A); } #elif LDBL_MANT_DIG <= 112 // quadruple precison @@ -374,9 +386,9 @@ void MatrixExponential<MatrixType>::computeUV(long double) pade13(m_M); } else { const long double maxnorm = 2.884233277829519311757165057717815L; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp<ScalingOp, MatrixType>(m_M, m_squarings); pade17(A); } #else |