From 1cd4279b03d5cb4a9ae16eef7af78b4af1003b8f Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 25 Aug 2012 01:09:20 +0800 Subject: Fix a lot in MatrixPower.h --- unsupported/test/matrix_power.cpp | 152 ++++++++++++++++++++++++++++++-------- 1 file changed, 122 insertions(+), 30 deletions(-) (limited to 'unsupported/test/matrix_power.cpp') diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index f3ef57157..80f65ebe4 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -23,7 +23,7 @@ void test2dRotation(double tol) 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"; + std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C, B) << '\n'; VERIFY(C.isApprox(B, T(tol))); } } @@ -43,44 +43,117 @@ void test2dHyperbolicRotation(double tol) B << ch, ish, -ish, ch; C = A.pow(angle); - std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C, B) << "\n"; + std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C, B) << '\n'; VERIFY(C.isApprox(B, T(tol))); } } template -void testExponentLaws(const MatrixType& m, double tol) +void testIntPowers(const MatrixType& m, double tol) { - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; + typedef typename MatrixType::RealScalar RealScalar; + const MatrixType m1 = MatrixType::Random(m.rows(), m.cols()); + const MatrixType identity = MatrixType::Identity(m.rows(), m.cols()); + const PartialPivLU solver(m1); + MatrixType m2, m3, m4; + + m3 = m1.pow(0); + m4 = m1.pow(0.); + std::cout << "testIntPower: i = 0 error powerm = " << relerr(identity, m3) << " " << relerr(identity, m4) << '\n'; + VERIFY(identity == m3 && identity == m4); + + m3 = m1.pow(1); + m4 = m1.pow(1.); + std::cout << "testIntPower: i = 1 error powerm = " << relerr(m1, m3) << " " << relerr(m1, m4) << '\n'; + VERIFY(m1 == m3 && m1 == m4); + + m2 = m1 * m1; + m3 = m1.pow(2); + m4 = m1.pow(2.); + std::cout << "testIntPower: i = 2 error powerm = " << relerr(m2, m3) << " " << relerr(m2, m4) << '\n'; + VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol))); + + for (int i = 3; i <= 20; i++) { + m2 *= m1; + m3 = m1.pow(i); + m4 = m1.pow(RealScalar(i)); + std::cout << "testIntPower: i = " << i << " error powerm = " << relerr(m2, m3) << " " << relerr (m2, m4) << '\n'; + VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol))); + } - typename MatrixType::Index rows = m.rows(); - typename MatrixType::Index cols = m.cols(); - MatrixType m1, m1x, m1y, m2, m3; - RealScalar x = internal::random(), y = internal::random(); - double err[3]; + m2 = solver.inverse(); + m3 = m1.pow(-1); + m4 = m1.pow(-1.); + std::cout << "testIntPower: i = -1 error powerm = " << relerr(m2, m3) << " " << relerr (m2, m4) << '\n'; + VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol))); - for(int i = 0; i < g_repeat; i++) { + for (int i = -2; i >= -20; i--) { + m2 = solver.solve(m2); + m3 = m1.pow(i); + m4 = m1.pow(RealScalar(i)); + std::cout << "testIntPower: i = " << i << " error powerm = " << relerr(m2, m3) << " " << relerr (m2, m4) << '\n'; + VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol))); + } +} + +template +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::run(m1, m.rows()); - m1x = m1.pow(x); - m1y = m1.pow(y); + x = internal::random(); + y = internal::random(); + m2 = m1.pow(x); + m3 = m1.pow(y); + + m4 = m1.pow(x + y); + m5 = m2 * m3; + std::cout << "testExponentLaws: error powerm = " << relerr(m4, m5); + VERIFY(m4.isApprox(m5, RealScalar(tol))); - m2 = m1.pow(x + y); - m3 = m1x * m1y; - err[0] = relerr(m2, m3); - VERIFY(m2.isApprox(m3, static_cast(tol))); + if (!NumTraits::IsComplex) { + m4 = m1.pow(x * y); + m5 = m2.pow(y); + std::cout << " " << relerr(m4, m5); + VERIFY(m4.isApprox(m5, RealScalar(tol))); + } - m2 = m1.pow(x * y); - m3 = m1x.pow(y); - err[1] = relerr(m2, m3); - VERIFY(m2.isApprox(m3, static_cast(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))); + } +} - 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(tol))); +template +void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double tol) +{ + typedef typename MatrixType::RealScalar RealScalar; + MatrixType m1; + VectorType v1, v2, v3; + RealScalar pReal; + signed char pInt; - std::cout << "testExponentLaws: error powerm = " << err[0] << " " << err[1] << " " << err[2] << "\n"; + for (int i = 0; i < g_repeat; i++) { + generateTestMatrix::run(m1, m.rows()); + v1 = VectorType::Random(v.rows(), v.cols()); + pReal = internal::random(); + pInt = rand(); + pInt >>= 2; + + v2 = m1.pow(pReal).eval() * v1; + v3 = m1.pow(pReal) * v1; + std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3); + VERIFY(v2.isApprox(v3, RealScalar(tol))); + + v2 = m1.pow(pInt).eval() * v1; + v3 = m1.pow(pInt) * v1; + std::cout << " " << relerr(v2, v3) << '\n'; + VERIFY(v2.isApprox(v3, RealScalar(tol)) || v2 == v3); } } @@ -88,17 +161,36 @@ void test_matrix_power() { CALL_SUBTEST_2(test2dRotation(1e-13)); CALL_SUBTEST_1(test2dRotation(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64 - CALL_SUBTEST_8(test2dRotation(1e-13)); + CALL_SUBTEST_9(test2dRotation(1e-13)); CALL_SUBTEST_2(test2dHyperbolicRotation(1e-14)); CALL_SUBTEST_1(test2dHyperbolicRotation(1e-5)); - CALL_SUBTEST_8(test2dHyperbolicRotation(1e-14)); + CALL_SUBTEST_9(test2dHyperbolicRotation(1e-14)); + + CALL_SUBTEST_2(testIntPowers(Matrix2d(), 1e-13)); + CALL_SUBTEST_7(testIntPowers(Matrix(), 1e-13)); + CALL_SUBTEST_3(testIntPowers(Matrix4cd(), 1e-13)); + CALL_SUBTEST_4(testIntPowers(MatrixXd(8,8), 1e-13)); + CALL_SUBTEST_1(testIntPowers(Matrix2f(), 1e-4)); + CALL_SUBTEST_5(testIntPowers(Matrix3cf(), 1e-4)); + CALL_SUBTEST_8(testIntPowers(Matrix4f(), 1e-4)); + CALL_SUBTEST_6(testIntPowers(MatrixXf(8,8), 1e-4)); + CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13)); CALL_SUBTEST_7(testExponentLaws(Matrix(), 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_8(testExponentLaws(Matrix4f(), 1e-4)); CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4)); - CALL_SUBTEST_9(testExponentLaws(Matrix(7,7), 1e-13)); + + CALL_SUBTEST_2(testMatrixVectorProduct(Matrix2d(), Vector2d(), 1e-13)); + CALL_SUBTEST_7(testMatrixVectorProduct(Matrix(), 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_10(testMatrixVectorProduct(Matrix(7,7), Matrix(), 1e-13)); } -- cgit v1.2.3