API change: ei_matrix_exponential(A) --> A.exp(), etc

As discussed on mailing list; see
http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2010/02/msg00190.html

Eigen/src/Core/MatrixBase.h

@@ -372,6 +372,16 @@ template<typename Derived> class MatrixBase
     template<typename OtherScalar>
     void applyOnTheRight(int p, int q, const PlanarRotation<OtherScalar>& j);
 
+///////// MatrixFunctions module /////////
+
+    typedef typename ei_stem_function<Scalar>::type StemFunction;
+    const MatrixExponentialReturnValue<Derived> exp() const;
+    const MatrixFunctionReturnValue<Derived> matrixFunction(StemFunction f) const;
+    const MatrixFunctionReturnValue<Derived> cosh() const;
+    const MatrixFunctionReturnValue<Derived> sinh() const;
+    const MatrixFunctionReturnValue<Derived> cos() const;
+    const MatrixFunctionReturnValue<Derived> sin() const;
+
 #ifdef EIGEN2_SUPPORT
     template<typename ProductDerived, typename Lhs, typename Rhs>
     Derived& operator+=(const Flagged<ProductBase<ProductDerived, Lhs,Rhs>, 0,

Eigen/src/Core/util/ForwardDeclarations.h

@@ -167,6 +167,17 @@ template<typename Scalar,int Dim> class Translation;
 template<typename Scalar> class UniformScaling;
 template<typename MatrixType,int Direction> class Homogeneous;
 
+// MatrixFunctions module
+template<typename Derived> struct MatrixExponentialReturnValue;
+template<typename Derived> struct MatrixFunctionReturnValue;
+template <typename Scalar>
+struct ei_stem_function
+{
+  typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
+  typedef ComplexScalar type(ComplexScalar, int);
+};
+
+
 #ifdef EIGEN2_SUPPORT
 template<typename ExpressionType> class Cwise;
 #endif

unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h

@@ -290,8 +290,8 @@ void MatrixExponential<MatrixType>::computeUV(double)
   * 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
+  * MatrixBase::exp() and most of the time this is the only way it is
-  * it is used.
+  * used.
   */
 template<typename Derived> struct MatrixExponentialReturnValue
 : public ReturnByValue<MatrixExponentialReturnValue<Derived> >
@@ -381,11 +381,10 @@ struct ei_traits<MatrixExponentialReturnValue<Derived> >
  * \c complex<float> or \c complex<double> .
  */
 template <typename Derived>
-MatrixExponentialReturnValue<Derived>
+const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
-ei_matrix_exponential(const MatrixBase<Derived> &M)
 {
-  ei_assert(M.rows() == M.cols());
+  ei_assert(rows() == cols());
-  return MatrixExponentialReturnValue<Derived>(M.derived());
+  return MatrixExponentialReturnValue<Derived>(derived());
 }
 
 #endif // EIGEN_MATRIX_EXPONENTIAL

unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h

@@ -59,7 +59,7 @@ class MatrixFunction
       * \param[out] result  the function \p f applied to \p A, as
       * specified in the constructor.
       *
-      * See ei_matrix_function() for details on how this computation
+      * See MatrixBase::matrixFunction() for details on how this computation
       * is implemented.
       */
     template <typename ResultType> 
@@ -486,8 +486,8 @@ typename MatrixFunction<MatrixType,1>::DynMatrixType MatrixFunction<MatrixType,1
   * This class holds the argument to the matrix function 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_function() and related functions and most of the time
+  * matrixBase::matrixFunction() and related functions and most of the
-  * this is the only way it is used.
+  * time this is the only way it is used.
   */
 template<typename Derived> class MatrixFunctionReturnValue
 : public ReturnByValue<MatrixFunctionReturnValue<Derived> >
@@ -533,6 +533,9 @@ struct ei_traits<MatrixFunctionReturnValue<Derived> >
 };
 
 
+/********** MatrixBase methods **********/
+
+
 /** \ingroup MatrixFunctions_Module
   *
   * \brief Compute a matrix function.
@@ -571,7 +574,7 @@ struct ei_traits<MatrixFunctionReturnValue<Derived> >
   *     \end{array} \right]. \f]
   * This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
   * the z-axis. This is the same example as used in the documentation
-  * of ei_matrix_exponential().
+  * of MatrixBase::exp().
   *
   * \include MatrixFunction.cpp
   * Output: \verbinclude MatrixFunction.out
@@ -580,16 +583,14 @@ struct ei_traits<MatrixFunctionReturnValue<Derived> >
   * \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);
+  * A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B);
   * \endcode
   */
 template <typename Derived>
-MatrixFunctionReturnValue<Derived>
+const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f) const
-ei_matrix_function(const MatrixBase<Derived> &M,
-		   typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f)
 {
-  ei_assert(M.rows() == M.cols());
+  ei_assert(rows() == cols());
-  return MatrixFunctionReturnValue<Derived>(M.derived(), f);
+  return MatrixFunctionReturnValue<Derived>(derived(), f);
 }
 
 /** \ingroup MatrixFunctions_Module
@@ -599,19 +600,17 @@ ei_matrix_function(const MatrixBase<Derived> &M,
   * \param[in]  M  a square matrix.
   * \returns  expression representing \f$ \sin(M) \f$.
   * 
-  * This function calls ei_matrix_function() with StdStemFunctions::sin().
+  * This function calls matrixFunction() with StdStemFunctions::sin().
   *
   * \include MatrixSine.cpp
   * Output: \verbinclude MatrixSine.out
   */
 template <typename Derived>
-MatrixFunctionReturnValue<Derived>
+const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
-ei_matrix_sin(const MatrixBase<Derived>& M)
 {
-  ei_assert(M.rows() == M.cols());
+  ei_assert(rows() == cols());
-  typedef typename ei_traits<Derived>::Scalar Scalar;
   typedef typename ei_stem_function<Scalar>::ComplexScalar ComplexScalar;
-  return MatrixFunctionReturnValue<Derived>(M.derived(), StdStemFunctions<ComplexScalar>::sin);
+  return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin);
 }
 
 /** \ingroup MatrixFunctions_Module
@@ -621,18 +620,16 @@ ei_matrix_sin(const MatrixBase<Derived>& M)
   * \param[in]  M  a square matrix.
   * \returns  expression representing \f$ \cos(M) \f$.
   * 
-  * This function calls ei_matrix_function() with StdStemFunctions::cos().
+  * This function calls matrixFunction() with StdStemFunctions::cos().
   *
   * \sa ei_matrix_sin() for an example.
   */
 template <typename Derived>
-MatrixFunctionReturnValue<Derived>
+const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
-ei_matrix_cos(const MatrixBase<Derived>& M)
 {
-  ei_assert(M.rows() == M.cols());
+  ei_assert(rows() == cols());
-  typedef typename ei_traits<Derived>::Scalar Scalar;
   typedef typename ei_stem_function<Scalar>::ComplexScalar ComplexScalar;
-  return MatrixFunctionReturnValue<Derived>(M.derived(), StdStemFunctions<ComplexScalar>::cos);
+  return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos);
 }
 
 /** \ingroup MatrixFunctions_Module
@@ -642,19 +639,17 @@ ei_matrix_cos(const MatrixBase<Derived>& M)
   * \param[in]  M  a square matrix.
   * \returns  expression representing \f$ \sinh(M) \f$
   * 
-  * This function calls ei_matrix_function() with StdStemFunctions::sinh().
+  * This function calls matrixFunction() with StdStemFunctions::sinh().
   *
   * \include MatrixSinh.cpp
   * Output: \verbinclude MatrixSinh.out
   */
 template <typename Derived>
-MatrixFunctionReturnValue<Derived>
+const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
-ei_matrix_sinh(const MatrixBase<Derived>& M)
 {
-  ei_assert(M.rows() == M.cols());
+  ei_assert(rows() == cols());
-  typedef typename ei_traits<Derived>::Scalar Scalar;
   typedef typename ei_stem_function<Scalar>::ComplexScalar ComplexScalar;
-  return MatrixFunctionReturnValue<Derived>(M.derived(), StdStemFunctions<ComplexScalar>::sinh);
+  return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh);
 }
 
 /** \ingroup MatrixFunctions_Module
@@ -664,18 +659,16 @@ ei_matrix_sinh(const MatrixBase<Derived>& M)
   * \param[in]  M  a square matrix.
   * \returns  expression representing \f$ \cosh(M) \f$
   * 
-  * This function calls ei_matrix_function() with StdStemFunctions::cosh().
+  * This function calls matrixFunction() with StdStemFunctions::cosh().
   *
   * \sa ei_matrix_sinh() for an example.
   */
 template <typename Derived>
-MatrixFunctionReturnValue<Derived>
+const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
-ei_matrix_cosh(const MatrixBase<Derived>& M)
 {
-  ei_assert(M.rows() == M.cols());
+  ei_assert(rows() == cols());
-  typedef typename ei_traits<Derived>::Scalar Scalar;
   typedef typename ei_stem_function<Scalar>::ComplexScalar ComplexScalar;
-  return MatrixFunctionReturnValue<Derived>(M.derived(), StdStemFunctions<ComplexScalar>::cosh);
+  return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cosh);
 }
 
 #endif // EIGEN_MATRIX_FUNCTION

unsupported/Eigen/src/MatrixFunctions/StemFunction.h

@@ -25,13 +25,6 @@
 #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.
   */

unsupported/doc/examples/MatrixExponential.cpp

@@ -12,7 +12,5 @@ int main()
        pi/4, 0,     0,
        0,    0,     0;
   std::cout << "The matrix A is:\n" << A << "\n\n";
-
+  std::cout << "The matrix exponential of A is:\n" << A.exp() << "\n\n";
-  MatrixXd B = ei_matrix_exponential(A);
-  std::cout << "The matrix exponential of A is:\n" << B << "\n\n";
 }

unsupported/doc/examples/MatrixFunction.cpp

@@ -19,5 +19,5 @@ int main()
 
   std::cout << "The matrix A is:\n" << A << "\n\n";
   std::cout << "The matrix exponential of A is:\n" 
-            << ei_matrix_function(A, expfn) << "\n\n";
+            << A.matrixFunction(expfn) << "\n\n";
 }

unsupported/doc/examples/MatrixSine.cpp

@@ -8,10 +8,10 @@ int main()
   MatrixXd A = MatrixXd::Random(3,3);
   std::cout << "A = \n" << A << "\n\n";
 
-  MatrixXd sinA = ei_matrix_sin(A);
+  MatrixXd sinA = A.sin();
   std::cout << "sin(A) = \n" << sinA << "\n\n";
 
-  MatrixXd cosA = ei_matrix_cos(A);
+  MatrixXd cosA = A.cos();
   std::cout << "cos(A) = \n" << cosA << "\n\n";
   
   // The matrix functions satisfy sin^2(A) + cos^2(A) = I, 

unsupported/doc/examples/MatrixSinh.cpp

@@ -8,10 +8,10 @@ int main()
   MatrixXf A = MatrixXf::Random(3,3);
   std::cout << "A = \n" << A << "\n\n";
 
-  MatrixXf sinhA = ei_matrix_sinh(A);
+  MatrixXf sinhA = A.sinh();
   std::cout << "sinh(A) = \n" << sinhA << "\n\n";
 
-  MatrixXf coshA = ei_matrix_cosh(A);
+  MatrixXf coshA = A.cosh();
   std::cout << "cosh(A) = \n" << coshA << "\n\n";
   
   // The matrix functions satisfy cosh^2(A) - sinh^2(A) = I, 

unsupported/test/matrix_exponential.cpp

@@ -57,11 +57,11 @@ void test2dRotation(double tol)
     angle = static_cast<T>(pow(10, i / 5. - 2));
     B << cos(angle), sin(angle), -sin(angle), cos(angle);
 
-    C = ei_matrix_function(angle*A, expfn);
+    C = (angle*A).matrixFunction(expfn);
     std::cout << "test2dRotation: i = " << i << "   error funm = " << relerr(C, B);
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
 
-    C = ei_matrix_exponential(angle*A);
+    C = (angle*A).exp();
     std::cout << "   error expm = " << relerr(C, B) << "\n";
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
   }
@@ -82,11 +82,11 @@ void test2dHyperbolicRotation(double tol)
     A << 0, angle*imagUnit, -angle*imagUnit, 0;
     B << ch, sh*imagUnit, -sh*imagUnit, ch;
 
-    C = ei_matrix_function(A, expfn);
+    C = A.matrixFunction(expfn);
     std::cout << "test2dHyperbolicRotation: i = " << i << "   error funm = " << relerr(C, B);
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
 
-    C = ei_matrix_exponential(A);
+    C = A.exp();
     std::cout << "   error expm = " << relerr(C, B) << "\n";
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
   }
@@ -106,11 +106,11 @@ void testPascal(double tol)
       for (int j=0; j<=i; j++)
     B(i,j) = static_cast<T>(binom(i,j));
 
-    C = ei_matrix_function(A, expfn);
+    C = A.matrixFunction(expfn);
     std::cout << "testPascal: size = " << size << "   error funm = " << relerr(C, B);
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
 
-    C = ei_matrix_exponential(A);
+    C = A.exp();
     std::cout << "   error expm = " << relerr(C, B) << "\n";
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
   }
@@ -132,11 +132,11 @@ void randomTest(const MatrixType& m, double tol)
   for(int i = 0; i < g_repeat; i++) {
     m1 = MatrixType::Random(rows, cols);
 
-    m2 = ei_matrix_function(m1, expfn) * ei_matrix_function(-m1, expfn);
+    m2 = m1.matrixFunction(expfn) * (-m1).matrixFunction(expfn);
     std::cout << "randomTest: error funm = " << relerr(identity, m2);
     VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol)));
 
-    m2 = ei_matrix_exponential(m1) * ei_matrix_exponential(-m1);
+    m2 = m1.exp() * (-m1).exp();
     std::cout << "   error expm = " << relerr(identity, m2) << "\n";
     VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol)));
   }

unsupported/test/matrix_function.cpp

@@ -114,8 +114,7 @@ void testMatrixExponential(const MatrixType& A)
   typedef typename NumTraits<Scalar>::Real RealScalar;
   typedef std::complex<RealScalar> ComplexScalar;
 
-  VERIFY_IS_APPROX(ei_matrix_exponential(A), 
+  VERIFY_IS_APPROX(A.exp(), A.matrixFunction(StdStemFunctions<ComplexScalar>::exp));
-		   ei_matrix_function(A, StdStemFunctions<ComplexScalar>::exp));
 }
 
 template<typename MatrixType>
@@ -123,10 +122,8 @@ void testHyperbolicFunctions(const MatrixType& A)
 {
   // Need to use absolute error because of possible cancellation when
   // adding/subtracting expA and expmA.
-  MatrixType expA = ei_matrix_exponential(A);
+  VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2);
-  MatrixType expmA = ei_matrix_exponential(-A);
+  VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2);
-  VERIFY_IS_APPROX_ABS(ei_matrix_sinh(A), (expA - expmA) / 2);
-  VERIFY_IS_APPROX_ABS(ei_matrix_cosh(A), (expA + expmA) / 2);
 }
 
 template<typename MatrixType>
@@ -143,15 +140,13 @@ void testGonioFunctions(const MatrixType& A)
 
   ComplexMatrix Ac = A.template cast<ComplexScalar>();
   
-  ComplexMatrix exp_iA = ei_matrix_exponential(imagUnit * Ac);
+  ComplexMatrix exp_iA = (imagUnit * Ac).exp();
-  ComplexMatrix exp_miA = ei_matrix_exponential(-imagUnit * Ac);
+  ComplexMatrix exp_miA = (-imagUnit * Ac).exp();
   
-  MatrixType sinA = ei_matrix_sin(A);
+  ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>();
-  ComplexMatrix sinAc = sinA.template cast<ComplexScalar>();
   VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
   
-  MatrixType cosA = ei_matrix_cos(A);
+  ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>();
-  ComplexMatrix cosAc = cosA.template cast<ComplexScalar>();
   VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2);
 }