From ba4e886376382ca29b5ae6b9f31b869805ab415a Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Tue, 28 Aug 2012 01:55:13 +0800 Subject: Tidy up and write dox. --- unsupported/test/matrix_power.cpp | 112 ++++++++++++-------------------------- 1 file changed, 35 insertions(+), 77 deletions(-) (limited to 'unsupported/test/matrix_power.cpp') diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 3c0e4f356..d891641f4 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -9,7 +9,7 @@ #include "matrix_functions.h" -template +template void test2dRotation(double tol) { Matrix A, B, C; @@ -28,7 +28,7 @@ void test2dRotation(double tol) } } -template +template void test2dHyperbolicRotation(double tol) { Matrix,2,2> A, B, C; @@ -48,55 +48,7 @@ void test2dHyperbolicRotation(double tol) } } -template -void testIntPowers(const MatrixType& m, double tol) -{ - 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.noalias() = 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))); - } - - 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 = -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 +template void testExponentLaws(const MatrixType& m, double tol) { typedef typename MatrixType::RealScalar RealScalar; @@ -129,31 +81,40 @@ void testExponentLaws(const MatrixType& m, double tol) } } -template +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; + RealScalar p; 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.noalias() = m1.pow(pReal).eval() * v1; - v3.noalias() = m1.pow(pReal) * v1; - std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3); - VERIFY(v2.isApprox(v3, RealScalar(tol))); - - v2.noalias() = m1.pow(pInt).eval() * v1; - v3.noalias() = m1.pow(pInt) * v1; - std::cout << " " << relerr(v2, v3) << '\n'; - VERIFY(v2.isApprox(v3, RealScalar(tol)) || v2 == v3); + p = internal::random(); + + 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 +void testAliasing(const MatrixType& m) +{ + typedef typename MatrixType::RealScalar RealScalar; + MatrixType m1, m2; + RealScalar p; + + for (int i = 0; i < g_repeat; i++) { + generateTestMatrix::run(m1, m.rows()); + p = internal::random(); + + m2 = m1.pow(p); + m1 = m1.pow(p); + VERIFY(m1 == m2); } } @@ -166,15 +127,6 @@ void test_matrix_power() CALL_SUBTEST_1(test2dHyperbolicRotation(1e-5)); 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)); @@ -192,5 +144,11 @@ void test_matrix_power() 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)); + CALL_SUBTEST_9(testMatrixVectorProduct(Matrix(7,7), Matrix(), 1e-13)); + + CALL_SUBTEST_7(testAliasing(Matrix())); + CALL_SUBTEST_3(testAliasing(Matrix4cd())); + CALL_SUBTEST_4(testAliasing(MatrixXd(8,8))); + CALL_SUBTEST_5(testAliasing(Matrix3cf())); + CALL_SUBTEST_6(testAliasing(MatrixXf(8,8))); } -- cgit v1.2.3