diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index bf2223a6e..91790c8d2 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -41,6 +41,15 @@ namespace Eigen { * \brief This module aims to provide various methods for the computation of * matrix functions. * + * %Matrix functions are defined as follows. Suppose that \f$ f \f$ + * is an entire function (that is, a function on the complex plane + * that is everywhere complex differentiable). Then its Taylor + * series + * \f[ f(0) + f'(0) x + \frac{f''(0)}{2} x^2 + \frac{f'''(0)}{3!} x^3 + \cdots \f] + * converges to \f$ f(x) \f$. In this case, we can define the matrix + * function by the same series: + * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f] + * * \code * #include * \endcode diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index b5f4e2b6f..fd1938a87 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -33,8 +33,8 @@ * * \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. + * \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] diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h index 43539f549..49326cd0e 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2009 Jitse Niesen +// Copyright (C) 2009, 2010 Jitse Niesen // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -25,12 +25,8 @@ #ifndef EIGEN_MATRIX_FUNCTION #define EIGEN_MATRIX_FUNCTION -template -struct ei_stem_function -{ - typedef std::complex::Real> ComplexScalar; - typedef ComplexScalar type(ComplexScalar, int); -}; +#include "StemFunction.h" +#include "MatrixFunctionAtomic.h" /** \ingroup MatrixFunctions_Module * @@ -43,14 +39,15 @@ struct ei_stem_function * This function computes \f$ f(A) \f$ and stores the result in the * matrix pointed to by \p result. * - * %Matrix functions are defined as follows. Suppose that \f$ f \f$ - * is an entire function (that is, a function on the complex plane - * that is everywhere complex differentiable). Then its Taylor - * series - * \f[ f(0) + f'(0) x + \frac{f''(0)}{2} x^2 + \frac{f'''(0)}{3!} x^3 + \cdots \f] - * converges to \f$ f(x) \f$. In this case, we can define the matrix - * function by the same series: - * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f] + * Suppose that \p M is a matrix whose entries have type \c Scalar. + * Then, the second argument, \p f, should be a function with prototype + * \code + * ComplexScalar f(ComplexScalar, int) + * \endcode + * where \c ComplexScalar = \c std::complex if \c Scalar is + * real (e.g., \c float or \c double) and \c ComplexScalar = + * \c Scalar if \c Scalar is complex. The return value of \c f(x,n) + * should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x. * * This routine uses the algorithm described in: * Philip Davies and Nicholas J. Higham, @@ -73,19 +70,21 @@ struct ei_stem_function * the z-axis. This is the same example as used in the documentation * of ei_matrix_exponential(). * - * Note that the function \c expfn is defined for complex numbers \c x, - * even though the matrix \c A is over the reals. - * * \include MatrixFunction.cpp * Output: \verbinclude MatrixFunction.out + * + * Note that the function \c expfn is defined for complex numbers + * \c x, even though the matrix \c A is over the reals. Instead of + * \c expfn, we could also have used StdStemFunctions::exp: + * \code + * ei_matrix_function(A, StdStemFunctions >::exp, &B); + * \endcode */ template EIGEN_STRONG_INLINE void ei_matrix_function(const MatrixBase& M, typename ei_stem_function::Scalar>::type f, typename MatrixBase::PlainMatrixType* result); -#include "MatrixFunctionAtomic.h" - /** \ingroup MatrixFunctions_Module * \brief Helper class for computing matrix functions. @@ -510,4 +509,94 @@ EIGEN_STRONG_INLINE void ei_matrix_function(const MatrixBase& M, MatrixFunction(M, f, result); } +/** \ingroup MatrixFunctions_Module + * + * \brief Compute the matrix sine. + * + * \param[in] M a square matrix. + * \param[out] result pointer to matrix in which to store the result, \f$ \sin(M) \f$ + * + * This function calls ei_matrix_function() with StdStemFunctions::sin(). + * + * \include MatrixSine.cpp + * Output: \verbinclude MatrixSine.out + */ +template +EIGEN_STRONG_INLINE void ei_matrix_sin(const MatrixBase& M, + typename MatrixBase::PlainMatrixType* result) +{ + ei_assert(M.rows() == M.cols()); + typedef typename MatrixBase::PlainMatrixType PlainMatrixType; + typedef typename ei_traits::Scalar Scalar; + typedef typename ei_stem_function::ComplexScalar ComplexScalar; + MatrixFunction(M, StdStemFunctions::sin, result); +} + +/** \ingroup MatrixFunctions_Module + * + * \brief Compute the matrix cosine. + * + * \param[in] M a square matrix. + * \param[out] result pointer to matrix in which to store the result, \f$ \cos(M) \f$ + * + * This function calls ei_matrix_function() with StdStemFunctions::cos(). + * + * \sa ei_matrix_sin() for an example. + */ +template +EIGEN_STRONG_INLINE void ei_matrix_cos(const MatrixBase& M, + typename MatrixBase::PlainMatrixType* result) +{ + ei_assert(M.rows() == M.cols()); + typedef typename MatrixBase::PlainMatrixType PlainMatrixType; + typedef typename ei_traits::Scalar Scalar; + typedef typename ei_stem_function::ComplexScalar ComplexScalar; + MatrixFunction(M, StdStemFunctions::cos, result); +} + +/** \ingroup MatrixFunctions_Module + * + * \brief Compute the matrix hyperbolic sine. + * + * \param[in] M a square matrix. + * \param[out] result pointer to matrix in which to store the result, \f$ \sinh(M) \f$ + * + * This function calls ei_matrix_function() with StdStemFunctions::sinh(). + * + * \include MatrixSinh.cpp + * Output: \verbinclude MatrixSinh.out + */ +template +EIGEN_STRONG_INLINE void ei_matrix_sinh(const MatrixBase& M, + typename MatrixBase::PlainMatrixType* result) +{ + ei_assert(M.rows() == M.cols()); + typedef typename MatrixBase::PlainMatrixType PlainMatrixType; + typedef typename ei_traits::Scalar Scalar; + typedef typename ei_stem_function::ComplexScalar ComplexScalar; + MatrixFunction(M, StdStemFunctions::sinh, result); +} + +/** \ingroup MatrixFunctions_Module + * + * \brief Compute the matrix hyberpolic cosine. + * + * \param[in] M a square matrix. + * \param[out] result pointer to matrix in which to store the result, \f$ \cosh(M) \f$ + * + * This function calls ei_matrix_function() with StdStemFunctions::cosh(). + * + * \sa ei_matrix_sinh() for an example. + */ +template +EIGEN_STRONG_INLINE void ei_matrix_cosh(const MatrixBase& M, + typename MatrixBase::PlainMatrixType* result) +{ + ei_assert(M.rows() == M.cols()); + typedef typename MatrixBase::PlainMatrixType PlainMatrixType; + typedef typename ei_traits::Scalar Scalar; + typedef typename ei_stem_function::ComplexScalar ComplexScalar; + MatrixFunction(M, StdStemFunctions::cosh, result); +} + #endif // EIGEN_MATRIX_FUNCTION diff --git a/unsupported/Eigen/src/MatrixFunctions/StemFunction.h b/unsupported/Eigen/src/MatrixFunctions/StemFunction.h new file mode 100644 index 000000000..90965c7dd --- /dev/null +++ b/unsupported/Eigen/src/MatrixFunctions/StemFunction.h @@ -0,0 +1,123 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Jitse Niesen +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_STEM_FUNCTION +#define EIGEN_STEM_FUNCTION + +template +struct ei_stem_function +{ + typedef std::complex::Real> ComplexScalar; + typedef ComplexScalar type(ComplexScalar, int); +}; + +/** \ingroup MatrixFunctions_Module + * \brief Stem functions corresponding to standard mathematical functions. + */ +template +class StdStemFunctions +{ + public: + + /** \brief The exponential function (and its derivatives). */ + static Scalar exp(Scalar x, int) + { + return std::exp(x); + } + + /** \brief Cosine (and its derivatives). */ + static Scalar cos(Scalar x, int n) + { + Scalar res; + switch (n % 4) { + case 0: + res = std::cos(x); + break; + case 1: + res = -std::sin(x); + break; + case 2: + res = -std::cos(x); + break; + case 3: + res = std::sin(x); + break; + } + return res; + } + + /** \brief Sine (and its derivatives). */ + static Scalar sin(Scalar x, int n) + { + Scalar res; + switch (n % 4) { + case 0: + res = std::sin(x); + break; + case 1: + res = std::cos(x); + break; + case 2: + res = -std::sin(x); + break; + case 3: + res = -std::cos(x); + break; + } + return res; + } + + /** \brief Hyperbolic cosine (and its derivatives). */ + static Scalar cosh(Scalar x, int n) + { + Scalar res; + switch (n % 2) { + case 0: + res = std::cosh(x); + break; + case 1: + res = std::sinh(x); + break; + } + return res; + } + + /** \brief Hyperbolic sine (and its derivatives). */ + static Scalar sinh(Scalar x, int n) + { + Scalar res; + switch (n % 2) { + case 0: + res = std::sinh(x); + break; + case 1: + res = std::cosh(x); + break; + } + return res; + } + +}; // end of class StdStemFunctions + +#endif // EIGEN_STEM_FUNCTION diff --git a/unsupported/doc/examples/MatrixSine.cpp b/unsupported/doc/examples/MatrixSine.cpp new file mode 100644 index 000000000..f8780ac92 --- /dev/null +++ b/unsupported/doc/examples/MatrixSine.cpp @@ -0,0 +1,21 @@ +#include + +using namespace Eigen; + +int main() +{ + MatrixXd A = MatrixXd::Random(3,3); + std::cout << "A = \n" << A << "\n\n"; + + MatrixXd sinA; + ei_matrix_sin(A, &sinA); + std::cout << "sin(A) = \n" << sinA << "\n\n"; + + MatrixXd cosA; + ei_matrix_cos(A, &cosA); + std::cout << "cos(A) = \n" << cosA << "\n\n"; + + // The matrix functions satisfy sin^2(A) + cos^2(A) = I, + // like the scalar functions. + std::cout << "sin^2(A) + cos^2(A) = \n" << sinA*sinA + cosA*cosA << "\n\n"; +} diff --git a/unsupported/doc/examples/MatrixSinh.cpp b/unsupported/doc/examples/MatrixSinh.cpp new file mode 100644 index 000000000..488d95652 --- /dev/null +++ b/unsupported/doc/examples/MatrixSinh.cpp @@ -0,0 +1,21 @@ +#include + +using namespace Eigen; + +int main() +{ + MatrixXf A = MatrixXf::Random(3,3); + std::cout << "A = \n" << A << "\n\n"; + + MatrixXf sinhA; + ei_matrix_sinh(A, &sinhA); + std::cout << "sinh(A) = \n" << sinhA << "\n\n"; + + MatrixXf coshA; + ei_matrix_cosh(A, &coshA); + std::cout << "cosh(A) = \n" << coshA << "\n\n"; + + // The matrix functions satisfy cosh^2(A) - sinh^2(A) = I, + // like the scalar functions. + std::cout << "cosh^2(A) - sinh^2(A) = \n" << coshA*coshA - sinhA*sinhA << "\n\n"; +} diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index 77758696d..339fba7cc 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -15,6 +15,7 @@ ei_add_test(NumericalDiff) ei_add_test(autodiff) ei_add_test(BVH) ei_add_test(matrix_exponential) +ei_add_test(matrix_function) ei_add_test(alignedvector3) ei_add_test(FFT) diff --git a/unsupported/test/matrix_function.cpp b/unsupported/test/matrix_function.cpp new file mode 100644 index 000000000..0eb30eecb --- /dev/null +++ b/unsupported/test/matrix_function.cpp @@ -0,0 +1,115 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Jitse Niesen +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#include "main.h" +#include + +template +void testMatrixExponential(const MatrixType& m) +{ + typedef typename ei_traits::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef std::complex ComplexScalar; + + const int rows = m.rows(); + const int cols = m.cols(); + + for (int i = 0; i < g_repeat; i++) { + MatrixType A = MatrixType::Random(rows, cols); + MatrixType expA1, expA2; + ei_matrix_exponential(A, &expA1); + ei_matrix_function(A, StdStemFunctions::exp, &expA2); + VERIFY_IS_APPROX(expA1, expA2); + } +} + +template +void testHyperbolicFunctions(const MatrixType& m) +{ + const int rows = m.rows(); + const int cols = m.cols(); + + for (int i = 0; i < g_repeat; i++) { + MatrixType A = MatrixType::Random(rows, cols); + MatrixType sinhA, coshA, expA; + ei_matrix_sinh(A, &sinhA); + ei_matrix_cosh(A, &coshA); + ei_matrix_exponential(A, &expA); + VERIFY_IS_APPROX(sinhA, (expA - expA.inverse())/2); + VERIFY_IS_APPROX(coshA, (expA + expA.inverse())/2); + } +} + +template +void testGonioFunctions(const MatrixType& m) +{ + typedef ei_traits Traits; + typedef typename Traits::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef std::complex ComplexScalar; + typedef Matrix ComplexMatrix; + + const int rows = m.rows(); + const int cols = m.cols(); + ComplexScalar imagUnit(0,1); + ComplexScalar two(2,0); + + for (int i = 0; i < g_repeat; i++) { + MatrixType A = MatrixType::Random(rows, cols); + ComplexMatrix Ac = A.template cast(); + + ComplexMatrix exp_iA; + ei_matrix_exponential(imagUnit * Ac, &exp_iA); + + MatrixType sinA; + ei_matrix_sin(A, &sinA); + ComplexMatrix sinAc = sinA.template cast(); + VERIFY_IS_APPROX(sinAc, (exp_iA - exp_iA.inverse()) / (two*imagUnit)); + + MatrixType cosA; + ei_matrix_cos(A, &cosA); + ComplexMatrix cosAc = cosA.template cast(); + VERIFY_IS_APPROX(cosAc, (exp_iA + exp_iA.inverse()) / 2); + } +} + +template +void testMatrixType(const MatrixType& m) +{ + testMatrixExponential(m); + testHyperbolicFunctions(m); + testGonioFunctions(m); +} + +void test_matrix_function() +{ + CALL_SUBTEST_1(testMatrixType(Matrix())); + CALL_SUBTEST_2(testMatrixType(Matrix3cf())); + CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8))); + CALL_SUBTEST_4(testMatrixType(Matrix2d())); + CALL_SUBTEST_5(testMatrixType(Matrix())); + CALL_SUBTEST_6(testMatrixType(Matrix4cd())); + CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13))); +}