1 files changed, 133 insertions, 0 deletions
diff --git a/eigen/unsupported/test/matrix_power.cpp b/eigen/unsupported/test/matrix_power.cpp
new file mode 100644
index 0000000..b9d513b
--- /dev/null
+++ b/eigen/unsupported/test/matrix_power.cpp
@@ -0,0 +1,133 @@
+// 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 MatrixType, int IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
+struct generateTriangularMatrix;
+
+// for real matrices, make sure none of the eigenvalues are negative
+template <typename MatrixType>
+struct generateTriangularMatrix<MatrixType,0>
+{
+  static void run(MatrixType& result, typename MatrixType::Index size)
+  {
+    result.resize(size, size);
+    result.template triangularView<Upper>() = MatrixType::Random(size, size);
+    for (typename MatrixType::Index i = 0; i < size; ++i)
+      result.coeffRef(i,i) = std::abs(result.coeff(i,i));
+  }
+};
+
+// for complex matrices, any matrix is fine
+template <typename MatrixType>
+struct generateTriangularMatrix<MatrixType,1>
+{
+  static void run(MatrixType& result, typename MatrixType::Index size)
+  {
+    result.resize(size, size);
+    result.template triangularView<Upper>() = MatrixType::Random(size, size);
+  }
+};
+
+template<typename T>
+void test2dRotation(double tol)
+{
+  Matrix<T,2,2> A, B, C;
+  T angle, c, s;
+
+  A << 0, 1, -1, 0;
+  MatrixPower<Matrix<T,2,2> > Apow(A);
+
+  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 = Apow(std::ldexp(angle,1) / M_PI);
+    std::cout << "test2dRotation: i = " << i << "   error powerm = " << relerr(C,B) << '\n';
+    VERIFY(C.isApprox(B, static_cast<T>(tol)));
+  }
+}
+
+template<typename T>
+void test2dHyperbolicRotation(double tol)
+{
+  Matrix<std::complex<T>,2,2> A, B, C;
+  T angle, ch = std::cosh((T)1);
+  std::complex<T> ish(0, std::sinh((T)1));
+
+  A << ch, ish, -ish, ch;
+  MatrixPower<Matrix<std::complex<T>,2,2> > Apow(A);
+
+  for (int i=0; i<=20; ++i) {
+    angle = std::ldexp(static_cast<T>(i-10), -1);
+    ch = std::cosh(angle);
+    ish = std::complex<T>(0, std::sinh(angle));
+    B << ch, ish, -ish, ch;
+
+    C = Apow(angle);
+    std::cout << "test2dHyperbolicRotation: i = " << i << "   error powerm = " << relerr(C,B) << '\n';
+    VERIFY(C.isApprox(B, static_cast<T>(tol)));
+  }
+}
+
+template<typename MatrixType>
+void testExponentLaws(const MatrixType& m, double tol)
+{
+  typedef typename MatrixType::RealScalar RealScalar;
+  MatrixType m1, m2, m3, m4, m5;
+  RealScalar x, y;
+
+  for (int i=0; i < g_repeat; ++i) {
+    generateTestMatrix<MatrixType>::run(m1, m.rows());
+    MatrixPower<MatrixType> mpow(m1);
+
+    x = internal::random<RealScalar>();
+    y = internal::random<RealScalar>();
+    m2 = mpow(x);
+    m3 = mpow(y);
+
+    m4 = mpow(x+y);
+    m5.noalias() = m2 * m3;
+    VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
+
+    m4 = mpow(x*y);
+    m5 = m2.pow(y);
+    VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
+
+    m4 = (std::abs(x) * m1).pow(y);
+    m5 = std::pow(std::abs(x), y) * m3;
+    VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
+  }
+}
+
+typedef Matrix<double,3,3,RowMajor>         Matrix3dRowMajor;
+typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
+ 
+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_9(test2dRotation<long double>(1e-13)); 
+  CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
+  CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
+  CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14));
+
+  CALL_SUBTEST_2(testExponentLaws(Matrix2d(),         1e-13));
+  CALL_SUBTEST_7(testExponentLaws(Matrix3dRowMajor(), 1e-13));
+  CALL_SUBTEST_3(testExponentLaws(Matrix4cd(),        1e-13));
+  CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8),      2e-12));
+  CALL_SUBTEST_1(testExponentLaws(Matrix2f(),         1e-4));
+  CALL_SUBTEST_5(testExponentLaws(Matrix3cf(),        1e-4));
+  CALL_SUBTEST_8(testExponentLaws(Matrix4f(),         1e-4));
+  CALL_SUBTEST_6(testExponentLaws(MatrixXf(2,2),      1e-3)); // see bug 614
+  CALL_SUBTEST_9(testExponentLaws(MatrixXe(7,7),      1e-13));
+}