Add support for matrix sine, cosine, sinh and cosh.

8 files changed, 401 insertions, 22 deletions

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 <unsupported/Eigen/MatrixFunctions>
   * \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 <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2009, 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // 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 <typename Scalar>
-struct ei_stem_function
-{
-  typedef std::complex<typename NumTraits<Scalar>::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<Scalar> 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<std::complex<double> >::exp, &B);
+  * \endcode
   */
 template <typename Derived>
 EIGEN_STRONG_INLINE void ei_matrix_function(const MatrixBase<Derived>& M, 
 					    typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f,
 					    typename MatrixBase<Derived>::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<Derived>& M,
   MatrixFunction<PlainMatrixType>(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 <typename Derived>
+EIGEN_STRONG_INLINE void ei_matrix_sin(const MatrixBase<Derived>& M, 
+				       typename MatrixBase<Derived>::PlainMatrixType* result)
+{
+  ei_assert(M.rows() == M.cols());
+  typedef typename MatrixBase<Derived>::PlainMatrixType PlainMatrixType;
+  typedef typename ei_traits<PlainMatrixType>::Scalar Scalar;
+  typedef typename ei_stem_function<Scalar>::ComplexScalar ComplexScalar;
+  MatrixFunction<PlainMatrixType>(M, StdStemFunctions<ComplexScalar>::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 <typename Derived>
+EIGEN_STRONG_INLINE void ei_matrix_cos(const MatrixBase<Derived>& M, 
+				       typename MatrixBase<Derived>::PlainMatrixType* result)
+{
+  ei_assert(M.rows() == M.cols());
+  typedef typename MatrixBase<Derived>::PlainMatrixType PlainMatrixType;
+  typedef typename ei_traits<PlainMatrixType>::Scalar Scalar;
+  typedef typename ei_stem_function<Scalar>::ComplexScalar ComplexScalar;
+  MatrixFunction<PlainMatrixType>(M, StdStemFunctions<ComplexScalar>::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 <typename Derived>
+EIGEN_STRONG_INLINE void ei_matrix_sinh(const MatrixBase<Derived>& M, 
+					typename MatrixBase<Derived>::PlainMatrixType* result)
+{
+  ei_assert(M.rows() == M.cols());
+  typedef typename MatrixBase<Derived>::PlainMatrixType PlainMatrixType;
+  typedef typename ei_traits<PlainMatrixType>::Scalar Scalar;
+  typedef typename ei_stem_function<Scalar>::ComplexScalar ComplexScalar;
+  MatrixFunction<PlainMatrixType>(M, StdStemFunctions<ComplexScalar>::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 <typename Derived>
+EIGEN_STRONG_INLINE void ei_matrix_cosh(const MatrixBase<Derived>& M, 
+					typename MatrixBase<Derived>::PlainMatrixType* result)
+{
+  ei_assert(M.rows() == M.cols());
+  typedef typename MatrixBase<Derived>::PlainMatrixType PlainMatrixType;
+  typedef typename ei_traits<PlainMatrixType>::Scalar Scalar;
+  typedef typename ei_stem_function<Scalar>::ComplexScalar ComplexScalar;
+  MatrixFunction<PlainMatrixType>(M, StdStemFunctions<ComplexScalar>::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 <jitse@maths.leeds.ac.uk>
+//
+// 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 <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_STEM_FUNCTION
+#define EIGEN_STEM_FUNCTION
+
+template <typename Scalar>
+struct ei_stem_function
+{
+  typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
+  typedef ComplexScalar type(ComplexScalar, int);
+};
+
+/** \ingroup MatrixFunctions_Module 
+  * \brief Stem functions corresponding to standard mathematical functions.
+  */
+template <typename Scalar>
+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 <unsupported/Eigen/MatrixFunctions>
+
+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 <unsupported/Eigen/MatrixFunctions>
+
+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 <jitse@maths.leeds.ac.uk>
+//
+// 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 <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+#include <unsupported/Eigen/MatrixFunctions>
+
+template<typename MatrixType>
+void testMatrixExponential(const MatrixType& m)
+{
+  typedef typename ei_traits<MatrixType>::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  typedef std::complex<RealScalar> 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<ComplexScalar>::exp, &expA2);
+    VERIFY_IS_APPROX(expA1, expA2);
+  }
+}
+
+template<typename MatrixType>
+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<typename MatrixType>
+void testGonioFunctions(const MatrixType& m)
+{
+  typedef ei_traits<MatrixType> Traits;
+  typedef typename Traits::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  typedef std::complex<RealScalar> ComplexScalar;
+  typedef Matrix<ComplexScalar, Traits::RowsAtCompileTime, 
+                 Traits::ColsAtCompileTime, MatrixType::Options> 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<ComplexScalar>();
+
+    ComplexMatrix exp_iA;
+    ei_matrix_exponential(imagUnit * Ac, &exp_iA);
+
+    MatrixType sinA;
+    ei_matrix_sin(A, &sinA);
+    ComplexMatrix sinAc = sinA.template cast<ComplexScalar>();
+    VERIFY_IS_APPROX(sinAc, (exp_iA - exp_iA.inverse()) / (two*imagUnit));
+
+    MatrixType cosA;
+    ei_matrix_cos(A, &cosA);
+    ComplexMatrix cosAc = cosA.template cast<ComplexScalar>();
+    VERIFY_IS_APPROX(cosAc, (exp_iA + exp_iA.inverse()) / 2);
+  }
+}
+
+template<typename MatrixType>
+void testMatrixType(const MatrixType& m)
+{
+  testMatrixExponential(m);
+  testHyperbolicFunctions(m);
+  testGonioFunctions(m);
+}
+
+void test_matrix_function()
+{
+  CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>()));
+  CALL_SUBTEST_2(testMatrixType(Matrix3cf()));
+  CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8)));
+  CALL_SUBTEST_4(testMatrixType(Matrix2d()));
+  CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>()));
+  CALL_SUBTEST_6(testMatrixType(Matrix4cd()));
+  CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13)));
+}