Add test for matrix power.

Use Christoph Hertzberg's suggestion to use exponent laws.

5 files changed, 155 insertions, 42 deletions

diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index b34b151b1..ff0137ec6 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -33,6 +33,7 @@ endif()
 
 ei_add_test(matrix_exponential)
 ei_add_test(matrix_function)
+ei_add_test(matrix_power)
 ei_add_test(matrix_square_root)
 ei_add_test(alignedvector3)
 ei_add_test(FFT)
diff --git a/unsupported/test/matrix_exponential.cpp b/unsupported/test/matrix_exponential.cpp
index 695472f91..50dec083d 100644
--- a/unsupported/test/matrix_exponential.cpp
+++ b/unsupported/test/matrix_exponential.cpp
@@ -7,8 +7,7 @@
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-#include "main.h"
-#include <unsupported/Eigen/MatrixFunctions>
+#include "matrix_functions.h"
 
 double binom(int n, int k)
 {
@@ -18,12 +17,6 @@ double binom(int n, int k)
   return res;
 }
 
-template <typename Derived, typename OtherDerived>
-double relerr(const MatrixBase<Derived>& A, const MatrixBase<OtherDerived>& B)
-{
-  return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum()));
-}
-
 template <typename T>
 T expfn(T x, int)
 {
@@ -109,8 +102,7 @@ void randomTest(const MatrixType& m, double tol)
   */
   typename MatrixType::Index rows = m.rows();
   typename MatrixType::Index cols = m.cols();
-  MatrixType m1(rows, cols), m2(rows, cols), m3(rows, cols),
-             identity = MatrixType::Identity(rows, rows);
+  MatrixType m1(rows, cols), m2(rows, cols), identity = MatrixType::Identity(rows, cols);
 
   typedef typename NumTraits<typename internal::traits<MatrixType>::Scalar>::Real RealScalar;
 
diff --git a/unsupported/test/matrix_functions.h b/unsupported/test/matrix_functions.h
new file mode 100644
index 000000000..5817caef6
--- /dev/null
+++ b/unsupported/test/matrix_functions.h
@@ -0,0 +1,47 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009-2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "main.h"
+#include <unsupported/Eigen/MatrixFunctions>
+
+template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
+struct generateTestMatrix;
+
+// for real matrices, make sure none of the eigenvalues are negative
+template <typename MatrixType>
+struct generateTestMatrix<MatrixType,0>
+{
+  static void run(MatrixType& result, typename MatrixType::Index size)
+  {
+    MatrixType mat = MatrixType::Random(size, size);
+    EigenSolver<MatrixType> es(mat);
+    typename EigenSolver<MatrixType>::EigenvalueType eivals = es.eigenvalues();
+    for (typename MatrixType::Index i = 0; i < size; ++i) {
+      if (eivals(i).imag() == 0 && eivals(i).real() < 0)
+	eivals(i) = -eivals(i);
+    }
+    result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real();
+  }
+};
+
+// for complex matrices, any matrix is fine
+template <typename MatrixType>
+struct generateTestMatrix<MatrixType,1>
+{
+  static void run(MatrixType& result, typename MatrixType::Index size)
+  {
+    result = MatrixType::Random(size, size);
+  }
+};
+
+template <typename Derived, typename OtherDerived>
+double relerr(const MatrixBase<Derived>& A, const MatrixBase<OtherDerived>& B)
+{
+  return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum()));
+}
diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp
new file mode 100644
index 000000000..f3ef57157
--- /dev/null
+++ b/unsupported/test/matrix_power.cpp
@@ -0,0 +1,104 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "matrix_functions.h"
+
+template <typename T>
+void test2dRotation(double tol)
+{
+  Matrix<T,2,2> A, B, C;
+  T angle, c, s;
+
+  A << 0, 1, -1, 0;
+  for (int i = 0; i <= 20; i++) {
+    angle = pow(10, (i-10) / 5.);
+    c = std::cos(angle);
+    s = std::sin(angle);
+    B << c, s, -s, c;
+
+    C = A.pow(std::ldexp(angle, 1) / M_PI);
+    std::cout << "test2dRotation: i = " << i << "   error powerm = " << relerr(C, B) << "\n";
+    VERIFY(C.isApprox(B, T(tol)));
+  }
+}
+
+template <typename T>
+void test2dHyperbolicRotation(double tol)
+{
+  Matrix<std::complex<T>,2,2> A, B, C;
+  T angle, ch = std::cosh(1);
+  std::complex<T> ish(0, std::sinh(1));
+
+  A << ch, ish, -ish, ch;
+  for (int i = 0; i <= 20; i++) {
+    angle = std::ldexp(T(i-10), -1);
+    ch = std::cosh(angle);
+    ish = std::complex<T>(0, std::sinh(angle));
+    B << ch, ish, -ish, ch;
+
+    C = A.pow(angle);
+    std::cout << "test2dHyperbolicRotation: i = " << i << "   error powerm = " << relerr(C, B) << "\n";
+    VERIFY(C.isApprox(B, T(tol)));
+  }
+}
+
+template <typename MatrixType>
+void testExponentLaws(const MatrixType& m, double tol)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+
+  typename MatrixType::Index rows = m.rows();
+  typename MatrixType::Index cols = m.cols();
+  MatrixType m1, m1x, m1y, m2, m3;
+  RealScalar x = internal::random<RealScalar>(), y = internal::random<RealScalar>();
+  double err[3];
+
+  for(int i = 0; i < g_repeat; i++) {
+    generateTestMatrix<MatrixType>::run(m1, m.rows());
+    m1x = m1.pow(x);
+    m1y = m1.pow(y);
+
+    m2 = m1.pow(x + y);
+    m3 = m1x * m1y;
+    err[0] = relerr(m2, m3);
+    VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol)));
+
+    m2 = m1.pow(x * y);
+    m3 = m1x.pow(y);
+    err[1] = relerr(m2, m3);
+    VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol)));
+
+    m2 = (std::abs(x) * m1).pow(y);
+    m3 = std::pow(std::abs(x), y) * m1y;
+    err[2] = relerr(m2, m3);
+    VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol)));
+
+    std::cout << "testExponentLaws: error powerm = " << err[0] << "  " << err[1] << "  " << err[2] << "\n";
+  }
+}
+
+void test_matrix_power()
+{
+  CALL_SUBTEST_2(test2dRotation<double>(1e-13));
+  CALL_SUBTEST_1(test2dRotation<float>(2e-5));  // was 1e-5, relaxed for clang 2.8 / linux / x86-64
+  CALL_SUBTEST_8(test2dRotation<long double>(1e-13)); 
+  CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
+  CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
+  CALL_SUBTEST_8(test2dHyperbolicRotation<long double>(1e-14));
+  CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13));
+  CALL_SUBTEST_7(testExponentLaws(Matrix<double,3,3,RowMajor>(), 1e-13));
+  CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13));
+  CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 1e-13));
+  CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4));
+  CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4));
+  CALL_SUBTEST_1(testExponentLaws(Matrix4f(), 1e-4));
+  CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4));
+  CALL_SUBTEST_9(testExponentLaws(Matrix<long double,Dynamic,Dynamic>(7,7), 1e-13));
+}
diff --git a/unsupported/test/matrix_square_root.cpp b/unsupported/test/matrix_square_root.cpp
index 508619a7a..ea541e1ea 100644
--- a/unsupported/test/matrix_square_root.cpp
+++ b/unsupported/test/matrix_square_root.cpp
@@ -7,38 +7,7 @@
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-#include "main.h"
-#include <unsupported/Eigen/MatrixFunctions>
-
-template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
-struct generateTestMatrix;
-
-// for real matrices, make sure none of the eigenvalues are negative
-template <typename MatrixType>
-struct generateTestMatrix<MatrixType,0>
-{
-  static void run(MatrixType& result, typename MatrixType::Index size)
-  {
-    MatrixType mat = MatrixType::Random(size, size);
-    EigenSolver<MatrixType> es(mat);
-    typename EigenSolver<MatrixType>::EigenvalueType eivals = es.eigenvalues();
-    for (typename MatrixType::Index i = 0; i < size; ++i) {
-      if (eivals(i).imag() == 0 && eivals(i).real() < 0)
-	eivals(i) = -eivals(i);
-    }
-    result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real();
-  }
-};
-
-// for complex matrices, any matrix is fine
-template <typename MatrixType>
-struct generateTestMatrix<MatrixType,1>
-{
-  static void run(MatrixType& result, typename MatrixType::Index size)
-  {
-    result = MatrixType::Random(size, size);
-  }
-};
+#include "matrix_functions.h"
 
 template<typename MatrixType>
 void testMatrixSqrt(const MatrixType& m)