From 25019f08366a027e783617b4b8bd86f00237a5ee Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Mon, 15 Feb 2010 19:17:25 +0000 Subject: [PATCH] Use ReturnByValue to return result of ei_matrix_exponential() . --- .../src/MatrixFunctions/MatrixExponential.h | 228 +++++++++++------- .../doc/examples/MatrixExponential.cpp | 3 +- unsupported/test/matrix_exponential.cpp | 13 +- unsupported/test/matrix_function.cpp | 14 +- 4 files changed, 159 insertions(+), 99 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index fd1938a87..147e21bc1 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -29,74 +29,32 @@ template Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); } #endif -/** \ingroup MatrixFunctions_Module - * - * \brief Compute the matrix exponential. - * - * \param[in] M matrix whose exponential is to be computed. - * \param[out] result pointer to the matrix in which to store the result. - * - * The matrix exponential of \f$ M \f$ is defined by - * \f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f] - * The matrix exponential can be used to solve linear ordinary - * differential equations: the solution of \f$ y' = My \f$ with the - * initial condition \f$ y(0) = y_0 \f$ is given by - * \f$ y(t) = \exp(M) y_0 \f$. - * - * The cost of the computation is approximately \f$ 20 n^3 \f$ for - * matrices of size \f$ n \f$. The number 20 depends weakly on the - * norm of the matrix. - * - * The matrix exponential is computed using the scaling-and-squaring - * method combined with Padé approximation. The matrix is first - * rescaled, then the exponential of the reduced matrix is computed - * approximant, and then the rescaling is undone by repeated - * squaring. The degree of the Padé approximant is chosen such - * that the approximation error is less than the round-off - * error. However, errors may accumulate during the squaring phase. - * - * Details of the algorithm can be found in: Nicholas J. Higham, "The - * scaling and squaring method for the matrix exponential revisited," - * SIAM J. %Matrix Anal. Applic., 26:1179–1193, - * 2005. - * - * Example: The following program checks that - * \f[ \exp \left[ \begin{array}{ccc} - * 0 & \frac14\pi & 0 \\ - * -\frac14\pi & 0 & 0 \\ - * 0 & 0 & 0 - * \end{array} \right] = \left[ \begin{array}{ccc} - * \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ - * \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ - * 0 & 0 & 1 - * \end{array} \right]. \f] - * This corresponds to a rotation of \f$ \frac14\pi \f$ radians around - * the z-axis. - * - * \include MatrixExponential.cpp - * Output: \verbinclude MatrixExponential.out - * - * \note \p M has to be a matrix of \c float, \c double, - * \c complex or \c complex . - */ -template -EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase &M, - typename MatrixBase::PlainMatrixType* result); /** \ingroup MatrixFunctions_Module * \brief Class for computing the matrix exponential. + * \tparam MatrixType type of the argument of the exponential, + * expected to be an instantiation of the Matrix class template. */ template class MatrixExponential { public: - /** \brief Compute the matrix exponential. - * - * \param M matrix whose exponential is to be computed. - * \param result pointer to the matrix in which to store the result. - */ - MatrixExponential(const MatrixType &M, MatrixType *result); + /** \brief Constructor. + * + * The class stores a reference to \p M, so it should not be + * changed (or destroyed) before compute() is called. + * + * \param[in] M matrix whose exponential is to be computed. + */ + MatrixExponential(const MatrixType &M); + + /** \brief Computes the matrix exponential. + * + * \param[out] result the matrix exponential of \p M in the constructor. + */ + template + void compute(ResultType &result); private: @@ -109,7 +67,7 @@ class MatrixExponential { * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. * - * \param A Argument of matrix exponential + * \param[in] A Argument of matrix exponential */ void pade3(const MatrixType &A); @@ -118,7 +76,7 @@ class MatrixExponential { * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. * - * \param A Argument of matrix exponential + * \param[in] A Argument of matrix exponential */ void pade5(const MatrixType &A); @@ -127,7 +85,7 @@ class MatrixExponential { * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. * - * \param A Argument of matrix exponential + * \param[in] A Argument of matrix exponential */ void pade7(const MatrixType &A); @@ -136,7 +94,7 @@ class MatrixExponential { * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. * - * \param A Argument of matrix exponential + * \param[in] A Argument of matrix exponential */ void pade9(const MatrixType &A); @@ -145,7 +103,7 @@ class MatrixExponential { * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. * - * \param A Argument of matrix exponential + * \param[in] A Argument of matrix exponential */ void pade13(const MatrixType &A); @@ -171,10 +129,10 @@ class MatrixExponential { void computeUV(float); typedef typename ei_traits::Scalar Scalar; - typedef typename NumTraits::Scalar>::Real RealScalar; + typedef typename NumTraits::Real RealScalar; - /** \brief Pointer to matrix whose exponential is to be computed. */ - const MatrixType* m_M; + /** \brief Reference to matrix whose exponential is to be computed. */ + const MatrixType& m_M; /** \brief Even-degree terms in numerator of Padé approximant. */ MatrixType m_U; @@ -199,8 +157,8 @@ class MatrixExponential { }; template -MatrixExponential::MatrixExponential(const MatrixType &M, MatrixType *result) : - m_M(&M), +MatrixExponential::MatrixExponential(const MatrixType &M) : + m_M(M), m_U(M.rows(),M.cols()), m_V(M.rows(),M.cols()), m_tmp1(M.rows(),M.cols()), @@ -208,13 +166,20 @@ MatrixExponential::MatrixExponential(const MatrixType &M, MatrixType m_Id(MatrixType::Identity(M.rows(), M.cols())), m_squarings(0), m_l1norm(static_cast(M.cwiseAbs().colwise().sum().maxCoeff())) +{ + /* empty body */ +} + +template +template +void MatrixExponential::compute(ResultType &result) { computeUV(RealScalar()); m_tmp1 = m_U + m_V; // numerator of Pade approximant m_tmp2 = -m_U + m_V; // denominator of Pade approximant - *result = m_tmp2.partialPivLu().solve(m_tmp1); + result = m_tmp2.partialPivLu().solve(m_tmp1); for (int i=0; i @@ -286,13 +251,13 @@ template void MatrixExponential::computeUV(float) { if (m_l1norm < 4.258730016922831e-001) { - pade3(*m_M); + pade3(m_M); } else if (m_l1norm < 1.880152677804762e+000) { - pade5(*m_M); + pade5(m_M); } else { const float maxnorm = 3.925724783138660f; m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm))); - MatrixType A = *m_M / std::pow(Scalar(2), Scalar(static_cast(m_squarings))); + MatrixType A = m_M / std::pow(Scalar(2), Scalar(static_cast(m_squarings))); pade7(A); } } @@ -301,27 +266,126 @@ template void MatrixExponential::computeUV(double) { if (m_l1norm < 1.495585217958292e-002) { - pade3(*m_M); + pade3(m_M); } else if (m_l1norm < 2.539398330063230e-001) { - pade5(*m_M); + pade5(m_M); } else if (m_l1norm < 9.504178996162932e-001) { - pade7(*m_M); + pade7(m_M); } else if (m_l1norm < 2.097847961257068e+000) { - pade9(*m_M); + pade9(m_M); } else { const double maxnorm = 5.371920351148152; m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm))); - MatrixType A = *m_M / std::pow(Scalar(2), Scalar(m_squarings)); + MatrixType A = m_M / std::pow(Scalar(2), Scalar(m_squarings)); pade13(A); } } +/** \ingroup MatrixFunctions_Module + * + * \brief Proxy for the matrix exponential of some matrix (expression). + * + * \tparam Derived Type of the argument to the matrix exponential. + * + * This class holds the argument to the matrix exponential until it + * is assigned or evaluated for some other reason (so the argument + * should not be changed in the meantime). It is the return type of + * ei_matrix_exponential() and most of the time this is the only way + * it is used. + */ +template struct MatrixExponentialReturnValue +: public ReturnByValue > +{ + public: + /** \brief Constructor. + * + * \param[in] src %Matrix (expression) forming the argument of the + * matrix exponential. + */ + MatrixExponentialReturnValue(const Derived& src) : m_src(src) { } + + /** \brief Compute the matrix exponential. + * + * \param[out] result the matrix exponential of \p src in the + * constructor. + */ + template + inline void evalTo(ResultType& result) const + { + const typename ei_eval::type srcEvaluated = m_src.eval(); + MatrixExponential me(srcEvaluated); + me.compute(result); + } + + int rows() const { return m_src.rows(); } + int cols() const { return m_src.cols(); } + + protected: + const Derived& m_src; +}; + +template +struct ei_traits > +{ + typedef typename Derived::PlainMatrixType ReturnMatrixType; +}; + +/** \ingroup MatrixFunctions_Module + * + * \brief Compute the matrix exponential. + * + * \param[in] M matrix whose exponential is to be computed. + * \returns expression representing the matrix exponential of \p M. + * + * The matrix exponential of \f$ M \f$ is defined by + * \f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f] + * The matrix exponential can be used to solve linear ordinary + * differential equations: the solution of \f$ y' = My \f$ with the + * initial condition \f$ y(0) = y_0 \f$ is given by + * \f$ y(t) = \exp(M) y_0 \f$. + * + * The cost of the computation is approximately \f$ 20 n^3 \f$ for + * matrices of size \f$ n \f$. The number 20 depends weakly on the + * norm of the matrix. + * + * The matrix exponential is computed using the scaling-and-squaring + * method combined with Padé approximation. The matrix is first + * rescaled, then the exponential of the reduced matrix is computed + * approximant, and then the rescaling is undone by repeated + * squaring. The degree of the Padé approximant is chosen such + * that the approximation error is less than the round-off + * error. However, errors may accumulate during the squaring phase. + * + * Details of the algorithm can be found in: Nicholas J. Higham, "The + * scaling and squaring method for the matrix exponential revisited," + * SIAM J. %Matrix Anal. Applic., 26:1179–1193, + * 2005. + * + * Example: The following program checks that + * \f[ \exp \left[ \begin{array}{ccc} + * 0 & \frac14\pi & 0 \\ + * -\frac14\pi & 0 & 0 \\ + * 0 & 0 & 0 + * \end{array} \right] = \left[ \begin{array}{ccc} + * \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ + * \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ + * 0 & 0 & 1 + * \end{array} \right]. \f] + * This corresponds to a rotation of \f$ \frac14\pi \f$ radians around + * the z-axis. + * + * \include MatrixExponential.cpp + * Output: \verbinclude MatrixExponential.out + * + * \note \p M has to be a matrix of \c float, \c double, + * \c complex or \c complex . + */ template -EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase &M, - typename MatrixBase::PlainMatrixType* result) +MatrixExponentialReturnValue +ei_matrix_exponential(const MatrixBase &M) { ei_assert(M.rows() == M.cols()); - MatrixExponential::PlainMatrixType>(M, result); + return MatrixExponentialReturnValue(M.derived()); } #endif // EIGEN_MATRIX_EXPONENTIAL diff --git a/unsupported/doc/examples/MatrixExponential.cpp b/unsupported/doc/examples/MatrixExponential.cpp index 7dc2396df..801ee4d01 100644 --- a/unsupported/doc/examples/MatrixExponential.cpp +++ b/unsupported/doc/examples/MatrixExponential.cpp @@ -12,7 +12,6 @@ int main() 0, 0, 0; std::cout << "The matrix A is:\n" << A << "\n\n"; - MatrixXd B; - ei_matrix_exponential(A, &B); + MatrixXd B = ei_matrix_exponential(A); std::cout << "The matrix exponential of A is:\n" << B << "\n\n"; } diff --git a/unsupported/test/matrix_exponential.cpp b/unsupported/test/matrix_exponential.cpp index a5b40adde..6150439c5 100644 --- a/unsupported/test/matrix_exponential.cpp +++ b/unsupported/test/matrix_exponential.cpp @@ -61,7 +61,7 @@ void test2dRotation(double tol) std::cout << "test2dRotation: i = " << i << " error funm = " << relerr(C, B); VERIFY(C.isApprox(B, static_cast(tol))); - ei_matrix_exponential(angle*A, &C); + C = ei_matrix_exponential(angle*A); std::cout << " error expm = " << relerr(C, B) << "\n"; VERIFY(C.isApprox(B, static_cast(tol))); } @@ -86,7 +86,7 @@ void test2dHyperbolicRotation(double tol) std::cout << "test2dHyperbolicRotation: i = " << i << " error funm = " << relerr(C, B); VERIFY(C.isApprox(B, static_cast(tol))); - ei_matrix_exponential(A, &C); + C = ei_matrix_exponential(A); std::cout << " error expm = " << relerr(C, B) << "\n"; VERIFY(C.isApprox(B, static_cast(tol))); } @@ -110,7 +110,7 @@ void testPascal(double tol) std::cout << "testPascal: size = " << size << " error funm = " << relerr(C, B); VERIFY(C.isApprox(B, static_cast(tol))); - ei_matrix_exponential(A, &C); + C = ei_matrix_exponential(A); std::cout << " error expm = " << relerr(C, B) << "\n"; VERIFY(C.isApprox(B, static_cast(tol))); } @@ -137,10 +137,9 @@ void randomTest(const MatrixType& m, double tol) std::cout << "randomTest: error funm = " << relerr(identity, m2 * m3); VERIFY(identity.isApprox(m2 * m3, static_cast(tol))); - ei_matrix_exponential(m1, &m2); - ei_matrix_exponential(-m1, &m3); - std::cout << " error expm = " << relerr(identity, m2 * m3) << "\n"; - VERIFY(identity.isApprox(m2 * m3, static_cast(tol))); + m2 = ei_matrix_exponential(m1) * ei_matrix_exponential(-m1); + std::cout << " error expm = " << relerr(identity, m2) << "\n"; + VERIFY(identity.isApprox(m2, static_cast(tol))); } } diff --git a/unsupported/test/matrix_function.cpp b/unsupported/test/matrix_function.cpp index de63937ad..25134f21d 100644 --- a/unsupported/test/matrix_function.cpp +++ b/unsupported/test/matrix_function.cpp @@ -100,10 +100,9 @@ void testMatrixExponential(const MatrixType& A) typedef std::complex ComplexScalar; for (int i = 0; i < g_repeat; i++) { - MatrixType expA1, expA2; - ei_matrix_exponential(A, &expA1); - ei_matrix_function(A, StdStemFunctions::exp, &expA2); - VERIFY_IS_APPROX(expA1, expA2); + MatrixType expA; + ei_matrix_function(A, StdStemFunctions::exp, &expA); + VERIFY_IS_APPROX(ei_matrix_exponential(A), expA); } } @@ -111,10 +110,10 @@ template void testHyperbolicFunctions(const MatrixType& A) { for (int i = 0; i < g_repeat; i++) { - MatrixType sinhA, coshA, expA; + MatrixType sinhA, coshA; ei_matrix_sinh(A, &sinhA); ei_matrix_cosh(A, &coshA); - ei_matrix_exponential(A, &expA); + MatrixType expA = ei_matrix_exponential(A); VERIFY_IS_APPROX(sinhA, (expA - expA.inverse())/2); VERIFY_IS_APPROX(coshA, (expA + expA.inverse())/2); } @@ -136,8 +135,7 @@ void testGonioFunctions(const MatrixType& A) for (int i = 0; i < g_repeat; i++) { ComplexMatrix Ac = A.template cast(); - ComplexMatrix exp_iA; - ei_matrix_exponential(imagUnit * Ac, &exp_iA); + ComplexMatrix exp_iA = ei_matrix_exponential(imagUnit * Ac); MatrixType sinA; ei_matrix_sin(A, &sinA);