Enable saving intermidiate (Schur decomposition) but disable unstable specialization for matrix power-matrix product.

1 files changed, 28 insertions, 74 deletions

diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp
index d891641f4..c5967d2eb 100644
--- a/unsupported/test/matrix_power.cpp
+++ b/unsupported/test/matrix_power.cpp
@@ -16,15 +16,17 @@ void test2dRotation(double tol)
   T angle, c, s;
 
   A << 0, 1, -1, 0;
-  for (int i = 0; i <= 20; i++) {
+  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 = 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)));
+    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)));
   }
 }
 
@@ -36,15 +38,17 @@ void test2dHyperbolicRotation(double tol)
   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);
+  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 = A.pow(angle);
-    std::cout << "test2dHyperbolicRotation: i = " << i << "   error powerm = " << relerr(C, B) << '\n';
-    VERIFY(C.isApprox(B, T(tol)));
+    C = Apow(angle);
+    std::cout << "test2dHyperbolicRotation: i = " << i << "   error powerm = " << relerr(C,B) << '\n';
+    VERIFY(C.isApprox(B, static_cast<T>(tol)));
   }
 }
 
@@ -55,71 +59,37 @@ void testExponentLaws(const MatrixType& m, double tol)
   MatrixType m1, m2, m3, m4, m5;
   RealScalar x, y;
 
-  for (int i = 0; i < g_repeat; i++) {
+  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 = m1.pow(x);
-    m3 = m1.pow(y);
+    m2 = mpow(x);
+    m3 = mpow(y);
 
-    m4 = m1.pow(x + y);
+    m4 = mpow(x+y);
     m5.noalias() = m2 * m3;
+
     std::cout << "testExponentLaws: error powerm = " << relerr(m4, m5);
-    VERIFY(m4.isApprox(m5, RealScalar(tol)));
+    VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
 
-    if (!NumTraits<typename MatrixType::Scalar>::IsComplex) {
-      m4 = m1.pow(x * y);
-      m5 = m2.pow(y);
-      std::cout << "   " << relerr(m4, m5);
-      VERIFY(m4.isApprox(m5, RealScalar(tol)));
-    }
+    m4 = mpow(x*y);
+    m5 = m2.pow(y);
+    std::cout << "   " << relerr(m4, m5);
+    VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
 
     m4 = (std::abs(x) * m1).pow(y);
     m5 = std::pow(std::abs(x), y) * m3;
     std::cout << "   " << relerr(m4, m5) << '\n';
-    VERIFY(m4.isApprox(m5, RealScalar(tol)));
-  }
-}
-
-template<typename MatrixType, typename VectorType>
-void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double tol)
-{
-  typedef typename MatrixType::RealScalar RealScalar;
-  MatrixType m1;
-  VectorType v1, v2, v3;
-  RealScalar p;
-
-  for (int i = 0; i < g_repeat; i++) {
-    generateTestMatrix<MatrixType>::run(m1, m.rows());
-    v1 = VectorType::Random(v.rows(), v.cols());
-    p = internal::random<RealScalar>();
-
-    v2.noalias() = m1.pow(p).eval() * v1;
-    v1 = m1.pow(p) * v1;
-    std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v1) << '\n';
-    VERIFY(v2.isApprox(v1, RealScalar(tol)));
-  }
-}
-
-template<typename MatrixType>
-void testAliasing(const MatrixType& m)
-{
-  typedef typename MatrixType::RealScalar RealScalar;
-  MatrixType m1, m2;
-  RealScalar p;
-
-  for (int i = 0; i < g_repeat; i++) {
-    generateTestMatrix<MatrixType>::run(m1, m.rows());
-    p = internal::random<RealScalar>();
-
-    m2 = m1.pow(p);
-    m1 = m1.pow(p);
-    VERIFY(m1 == m2);
+    VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
   }
 }
 
 void test_matrix_power()
 {
+  typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
+
   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)); 
@@ -135,20 +105,4 @@ void test_matrix_power()
   CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4));
   CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4));
   CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4));
-
-  CALL_SUBTEST_2(testMatrixVectorProduct(Matrix2d(), Vector2d(), 1e-13));
-  CALL_SUBTEST_7(testMatrixVectorProduct(Matrix<double,3,3,RowMajor>(), Vector3d(), 1e-13));
-  CALL_SUBTEST_3(testMatrixVectorProduct(Matrix4cd(), Vector4cd(), 1e-13));
-  CALL_SUBTEST_4(testMatrixVectorProduct(MatrixXd(8,8), MatrixXd(8,2), 1e-13));
-  CALL_SUBTEST_1(testMatrixVectorProduct(Matrix2f(), Vector2f(), 1e-4));
-  CALL_SUBTEST_5(testMatrixVectorProduct(Matrix3cf(), Vector3cf(), 1e-4));
-  CALL_SUBTEST_8(testMatrixVectorProduct(Matrix4f(), Vector4f(), 1e-4));
-  CALL_SUBTEST_6(testMatrixVectorProduct(MatrixXf(8,8), VectorXf(8), 1e-4));
-  CALL_SUBTEST_9(testMatrixVectorProduct(Matrix<long double,Dynamic,Dynamic>(7,7), Matrix<long double,7,9>(), 1e-13));
-
-  CALL_SUBTEST_7(testAliasing(Matrix<double,3,3,RowMajor>()));
-  CALL_SUBTEST_3(testAliasing(Matrix4cd()));
-  CALL_SUBTEST_4(testAliasing(MatrixXd(8,8)));
-  CALL_SUBTEST_5(testAliasing(Matrix3cf()));
-  CALL_SUBTEST_6(testAliasing(MatrixXf(8,8)));
 }