diff options
Diffstat (limited to 'unsupported/test/matrix_power.cpp')
-rw-r--r-- | unsupported/test/matrix_power.cpp | 116 |
1 files changed, 49 insertions, 67 deletions
diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 3a3f01464..95c63c574 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((T)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,64 @@ 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); + 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))); + VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol))); } } template<typename MatrixType, typename VectorType> -void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double tol) +void testProduct(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++) { + for (int i=0; i<g_repeat; ++i) { generateTestMatrix<MatrixType>::run(m1, m.rows()); + MatrixPower<MatrixType> mpow(m1); + 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))); + v2.noalias() = mpow(p) * v1; + v3.noalias() = mpow(p).eval() * v1; + std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3) << '\n'; + VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol))); } } -template<typename MatrixType> -void testAliasing(const MatrixType& m) +template<typename MatrixType, typename VectorType> +void testMatrixVector(const MatrixType& m, const VectorType& v, double tol) { - 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); - } + testExponentLaws(m,tol); + testProduct(m,v,tol); } void test_matrix_power() { + typedef Matrix<double,3,3,RowMajor> Matrix3dRowMajor; + typedef Matrix<long double,Dynamic,Dynamic> MatrixXe; + typedef Matrix<long double,Dynamic,1> VectorXe; + 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)); @@ -127,28 +124,13 @@ void test_matrix_power() 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(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_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))); + CALL_SUBTEST_2(testMatrixVector(Matrix2d(), Vector2d(), 1e-13)); + CALL_SUBTEST_7(testMatrixVector(Matrix3dRowMajor(), MatrixXd(3,5), 1e-13)); + CALL_SUBTEST_3(testMatrixVector(Matrix4cd(), Vector4cd(), 1e-13)); + CALL_SUBTEST_4(testMatrixVector(MatrixXd(8,8), VectorXd(8), 1e-13)); + CALL_SUBTEST_1(testMatrixVector(Matrix2f(), Vector2f(), 1e-4)); + CALL_SUBTEST_5(testMatrixVector(Matrix3cf(), Vector3cf(), 1e-4)); + CALL_SUBTEST_8(testMatrixVector(Matrix4f(), Vector4f(), 1e-4)); + CALL_SUBTEST_6(testMatrixVector(MatrixXf(8,8), VectorXf(8), 1e-4)); + CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7), VectorXe(7), 1e-13)); } |