From 483eb2f56657e8e7f419ab1a4fab8dce9ade8609 Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sat, 27 Apr 2024 20:24:20 +0200 Subject: Adding upstream version 14.2.21. Signed-off-by: Daniel Baumann --- src/boost/libs/multiprecision/test/eigen.hpp | 627 +++++++++++++++++++++++++++ 1 file changed, 627 insertions(+) create mode 100644 src/boost/libs/multiprecision/test/eigen.hpp (limited to 'src/boost/libs/multiprecision/test/eigen.hpp') diff --git a/src/boost/libs/multiprecision/test/eigen.hpp b/src/boost/libs/multiprecision/test/eigen.hpp new file mode 100644 index 00000000..c01b949a --- /dev/null +++ b/src/boost/libs/multiprecision/test/eigen.hpp @@ -0,0 +1,627 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include "test.hpp" + +using namespace Eigen; + +template +struct related_number +{ + typedef T type; +}; + +template +void example1() +{ + // expected results first: + Matrix r1, r2; + r1 << 3, 5, 4, 8; + r2 << -1, -1, 2, 0; + Matrix r3; + r3 << -1, -4, -6; + + Matrix a; + a << 1, 2, 3, 4; + Matrix b(2, 2); + b << 2, 3, 1, 4; + std::cout << "a + b =\n" + << a + b << std::endl; + BOOST_CHECK_EQUAL(a + b, r1); + std::cout << "a - b =\n" + << a - b << std::endl; + BOOST_CHECK_EQUAL(a - b, r2); + std::cout << "Doing a += b;" << std::endl; + a += b; + std::cout << "Now a =\n" + << a << std::endl; + Matrix v(1, 2, 3); + Matrix w(1, 0, 0); + std::cout << "-v + w - v =\n" + << -v + w - v << std::endl; + BOOST_CHECK_EQUAL(-v + w - v, r3); +} + +template +void example2() +{ + Matrix a; + a << 1, 2, 3, 4; + Matrix v(1, 2, 3); + std::cout << "a * 2.5 =\n" + << a * 2.5 << std::endl; + std::cout << "0.1 * v =\n" + << 0.1 * v << std::endl; + std::cout << "Doing v *= 2;" << std::endl; + v *= 2; + std::cout << "Now v =\n" + << v << std::endl; + Num n(4); + std::cout << "Doing v *= Num;" << std::endl; + v *= n; + std::cout << "Now v =\n" + << v << std::endl; + typedef typename related_number::type related_type; + related_type r(6); + std::cout << "Doing v *= RelatedType;" << std::endl; + v *= r; + std::cout << "Now v =\n" + << v << std::endl; + std::cout << "RelatedType * v =\n" + << r * v << std::endl; + std::cout << "Doing v *= RelatedType^2;" << std::endl; + v *= r * r; + std::cout << "Now v =\n" + << v << std::endl; + std::cout << "RelatedType^2 * v =\n" + << r * r * v << std::endl; +} + +template +void example3() +{ + using namespace std; + Matrix a = Matrix::Random(2, 2); + cout << "Here is the matrix a\n" + << a << endl; + cout << "Here is the matrix a^T\n" + << a.transpose() << endl; + cout << "Here is the conjugate of a\n" + << a.conjugate() << endl; + cout << "Here is the matrix a^*\n" + << a.adjoint() << endl; +} + +template +void example4() +{ + Matrix mat; + mat << 1, 2, + 3, 4; + Matrix u(-1, 1), v(2, 0); + std::cout << "Here is mat*mat:\n" + << mat * mat << std::endl; + std::cout << "Here is mat*u:\n" + << mat * u << std::endl; + std::cout << "Here is u^T*mat:\n" + << u.transpose() * mat << std::endl; + std::cout << "Here is u^T*v:\n" + << u.transpose() * v << std::endl; + std::cout << "Here is u*v^T:\n" + << u * v.transpose() << std::endl; + std::cout << "Let's multiply mat by itself" << std::endl; + mat = mat * mat; + std::cout << "Now mat is mat:\n" + << mat << std::endl; +} + +template +void example5() +{ + using namespace std; + Matrix v(1, 2, 3); + Matrix w(0, 1, 2); + cout << "Dot product: " << v.dot(w) << endl; + Num dp = v.adjoint() * w; // automatic conversion of the inner product to a scalar + cout << "Dot product via a matrix product: " << dp << endl; + cout << "Cross product:\n" + << v.cross(w) << endl; +} + +template +void example6() +{ + using namespace std; + Matrix mat; + mat << 1, 2, + 3, 4; + cout << "Here is mat.sum(): " << mat.sum() << endl; + cout << "Here is mat.prod(): " << mat.prod() << endl; + cout << "Here is mat.mean(): " << mat.mean() << endl; + cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl; + cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl; + cout << "Here is mat.trace(): " << mat.trace() << endl; +} + +template +void example7() +{ + using namespace std; + + Array m(2, 2); + + // assign some values coefficient by coefficient + m(0, 0) = 1.0; + m(0, 1) = 2.0; + m(1, 0) = 3.0; + m(1, 1) = m(0, 1) + m(1, 0); + + // print values to standard output + cout << m << endl + << endl; + + // using the comma-initializer is also allowed + m << 1.0, 2.0, + 3.0, 4.0; + + // print values to standard output + cout << m << endl; +} + +template +void example8() +{ + using namespace std; + Array a(3, 3); + Array b(3, 3); + a << 1, 2, 3, + 4, 5, 6, + 7, 8, 9; + b << 1, 2, 3, + 1, 2, 3, + 1, 2, 3; + + // Adding two arrays + cout << "a + b = " << endl + << a + b << endl + << endl; + // Subtracting a scalar from an array + cout << "a - 2 = " << endl + << a - 2 << endl; +} + +template +void example9() +{ + using namespace std; + Array a(2, 2); + Array b(2, 2); + a << 1, 2, + 3, 4; + b << 5, 6, + 7, 8; + cout << "a * b = " << endl + << a * b << endl; +} + +template +void example10() +{ + using namespace std; + Array a = Array::Random(5); + a *= 2; + cout << "a =" << endl + << a << endl; + cout << "a.abs() =" << endl + << a.abs() << endl; + cout << "a.abs().sqrt() =" << endl + << a.abs().sqrt() << endl; + cout << "a.min(a.abs().sqrt()) =" << endl + << (a.min)(a.abs().sqrt()) << endl; +} + +template +void example11() +{ + using namespace std; + Matrix m(2, 2); + Matrix n(2, 2); + Matrix result(2, 2); + m << 1, 2, + 3, 4; + n << 5, 6, + 7, 8; + result = m * n; + cout << "-- Matrix m*n: --" << endl + << result << endl + << endl; + result = m.array() * n.array(); + cout << "-- Array m*n: --" << endl + << result << endl + << endl; + result = m.cwiseProduct(n); + cout << "-- With cwiseProduct: --" << endl + << result << endl + << endl; + result = m.array() + 4; + cout << "-- Array m + 4: --" << endl + << result << endl + << endl; +} + +template +void example12() +{ + using namespace std; + Matrix m(2, 2); + Matrix n(2, 2); + Matrix result(2, 2); + m << 1, 2, + 3, 4; + n << 5, 6, + 7, 8; + + result = (m.array() + 4).matrix() * m; + cout << "-- Combination 1: --" << endl + << result << endl + << endl; + result = (m.array() * n.array()).matrix() * m; + cout << "-- Combination 2: --" << endl + << result << endl + << endl; +} + +template +void example13() +{ + using namespace std; + Matrix m(4, 4); + m << 1, 2, 3, 4, + 5, 6, 7, 8, + 9, 10, 11, 12, + 13, 14, 15, 16; + cout << "Block in the middle" << endl; + cout << m.template block<2, 2>(1, 1) << endl + << endl; + for (int i = 1; i <= 3; ++i) + { + cout << "Block of size " << i << "x" << i << endl; + cout << m.block(0, 0, i, i) << endl + << endl; + } +} + +template +void example14() +{ + using namespace std; + Array m; + m << 1, 2, + 3, 4; + Array a = Array::Constant(0.6); + cout << "Here is the array a:" << endl + << a << endl + << endl; + a.template block<2, 2>(1, 1) = m; + cout << "Here is now a with m copied into its central 2x2 block:" << endl + << a << endl + << endl; + a.block(0, 0, 2, 3) = a.block(2, 1, 2, 3); + cout << "Here is now a with bottom-right 2x3 block copied into top-left 2x2 block:" << endl + << a << endl + << endl; +} + +template +void example15() +{ + using namespace std; + Eigen::Matrix m(3, 3); + m << 1, 2, 3, + 4, 5, 6, + 7, 8, 9; + cout << "Here is the matrix m:" << endl + << m << endl; + cout << "2nd Row: " << m.row(1) << endl; + m.col(2) += 3 * m.col(0); + cout << "After adding 3 times the first column into the third column, the matrix m is:\n"; + cout << m << endl; +} + +template +void example16() +{ + using namespace std; + Matrix m; + m << 1, 2, 3, 4, + 5, 6, 7, 8, + 9, 10, 11, 12, + 13, 14, 15, 16; + cout << "m.leftCols(2) =" << endl + << m.leftCols(2) << endl + << endl; + cout << "m.bottomRows<2>() =" << endl + << m.template bottomRows<2>() << endl + << endl; + m.topLeftCorner(1, 3) = m.bottomRightCorner(3, 1).transpose(); + cout << "After assignment, m = " << endl + << m << endl; +} + +template +void example17() +{ + using namespace std; + Array v(6); + v << 1, 2, 3, 4, 5, 6; + cout << "v.head(3) =" << endl + << v.head(3) << endl + << endl; + cout << "v.tail<3>() = " << endl + << v.template tail<3>() << endl + << endl; + v.segment(1, 4) *= 2; + cout << "after 'v.segment(1,4) *= 2', v =" << endl + << v << endl; +} + +template +void example18() +{ + using namespace std; + Matrix mat; + mat << 1, 2, + 3, 4; + cout << "Here is mat.sum(): " << mat.sum() << endl; + cout << "Here is mat.prod(): " << mat.prod() << endl; + cout << "Here is mat.mean(): " << mat.mean() << endl; + cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl; + cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl; + cout << "Here is mat.trace(): " << mat.trace() << endl; + + BOOST_CHECK_EQUAL(mat.sum(), 10); + BOOST_CHECK_EQUAL(mat.prod(), 24); + BOOST_CHECK_EQUAL(mat.mean(), Num(5) / 2); + BOOST_CHECK_EQUAL(mat.minCoeff(), 1); + BOOST_CHECK_EQUAL(mat.maxCoeff(), 4); + BOOST_CHECK_EQUAL(mat.trace(), 5); +} + +template +void example18a() +{ + using namespace std; + Matrix mat; + mat << 1, 2, + 3, 4; + cout << "Here is mat.sum(): " << mat.sum() << endl; + cout << "Here is mat.prod(): " << mat.prod() << endl; + cout << "Here is mat.mean(): " << mat.mean() << endl; + //cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl; + //cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl; + cout << "Here is mat.trace(): " << mat.trace() << endl; +} + +template +void example19() +{ + using namespace std; + Matrix v(2); + Matrix m(2, 2), n(2, 2); + + v << -1, + 2; + + m << 1, -2, + -3, 4; + cout << "v.squaredNorm() = " << v.squaredNorm() << endl; + cout << "v.norm() = " << v.norm() << endl; + cout << "v.lpNorm<1>() = " << v.template lpNorm<1>() << endl; + cout << "v.lpNorm() = " << v.template lpNorm() << endl; + cout << endl; + cout << "m.squaredNorm() = " << m.squaredNorm() << endl; + cout << "m.norm() = " << m.norm() << endl; + cout << "m.lpNorm<1>() = " << m.template lpNorm<1>() << endl; + cout << "m.lpNorm() = " << m.template lpNorm() << endl; +} + +template +void example20() +{ + using namespace std; + Matrix A; + Matrix b; + A << 1, 2, 3, 4, 5, 6, 7, 8, 10; + b << 3, 3, 4; + cout << "Here is the matrix A:\n" + << A << endl; + cout << "Here is the vector b:\n" + << b << endl; + Matrix x = A.colPivHouseholderQr().solve(b); + cout << "The solution is:\n" + << x << endl; +} + +template +void example21() +{ + using namespace std; + Matrix A, b; + A << 2, -1, -1, 3; + b << 1, 2, 3, 1; + cout << "Here is the matrix A:\n" + << A << endl; + cout << "Here is the right hand side b:\n" + << b << endl; + Matrix x = A.ldlt().solve(b); + cout << "The solution is:\n" + << x << endl; +} + +template +void example22() +{ + using namespace std; + Matrix A = Matrix::Random(100, 100); + Matrix b = Matrix::Random(100, 50); + Matrix x = A.fullPivLu().solve(b); + Matrix axmb = A * x - b; + double relative_error = static_cast(abs(axmb.norm() / b.norm())); // norm() is L2 norm + cout << "norm1 = " << axmb.norm() << endl; + cout << "norm2 = " << b.norm() << endl; + cout << "The relative error is:\n" + << relative_error << endl; +} + +template +void example23() +{ + using namespace std; + Matrix A; + A << 1, 2, 2, 3; + cout << "Here is the matrix A:\n" + << A << endl; + SelfAdjointEigenSolver > eigensolver(A); + if (eigensolver.info() != Success) + { + std::cout << "Eigenvalue solver failed!" << endl; + } + else + { + cout << "The eigenvalues of A are:\n" + << eigensolver.eigenvalues() << endl; + cout << "Here's a matrix whose columns are eigenvectors of A \n" + << "corresponding to these eigenvalues:\n" + << eigensolver.eigenvectors() << endl; + } +} + +template +void example24() +{ + using namespace std; + Matrix A; + A << 1, 2, 1, + 2, 1, 0, + -1, 1, 2; + cout << "Here is the matrix A:\n" + << A << endl; + cout << "The determinant of A is " << A.determinant() << endl; + cout << "The inverse of A is:\n" + << A.inverse() << endl; +} + +template +void test_integer_type() +{ + example1(); + //example2(); + example18(); +} + +template +void test_float_type() +{ + std::cout << "Epsilon = " << Eigen::NumTraits::epsilon() << std::endl; + std::cout << "Dummy Prec = " << Eigen::NumTraits::dummy_precision() << std::endl; + std::cout << "Highest = " << Eigen::NumTraits::highest() << std::endl; + std::cout << "Lowest = " << Eigen::NumTraits::lowest() << std::endl; + std::cout << "Digits10 = " << Eigen::NumTraits::digits10() << std::endl; + + example1(); + example2(); + example4(); + example5(); + example6(); + example7(); + example8(); + example9(); + example10(); + example11(); + example12(); + example13(); + example14(); + example15(); + example16(); + example17(); + /* + example18(); + example19(); + example20(); + example21(); + example22(); + example23(); + example24(); + */ +} + +template +void test_float_type_2() +{ + std::cout << "Epsilon = " << Eigen::NumTraits::epsilon() << std::endl; + std::cout << "Dummy Prec = " << Eigen::NumTraits::dummy_precision() << std::endl; + std::cout << "Highest = " << Eigen::NumTraits::highest() << std::endl; + std::cout << "Lowest = " << Eigen::NumTraits::lowest() << std::endl; + std::cout << "Digits10 = " << Eigen::NumTraits::digits10() << std::endl; + + example18(); + example19(); + example20(); + example21(); + + //example22(); + //example23(); + //example24(); +} + +template +void test_float_type_3() +{ + std::cout << "Epsilon = " << Eigen::NumTraits::epsilon() << std::endl; + std::cout << "Dummy Prec = " << Eigen::NumTraits::dummy_precision() << std::endl; + std::cout << "Highest = " << Eigen::NumTraits::highest() << std::endl; + std::cout << "Lowest = " << Eigen::NumTraits::lowest() << std::endl; + std::cout << "Digits10 = " << Eigen::NumTraits::digits10() << std::endl; + + example22(); + example23(); + example24(); +} + +template +void test_complex_type() +{ + std::cout << "Epsilon = " << Eigen::NumTraits::epsilon() << std::endl; + std::cout << "Dummy Prec = " << Eigen::NumTraits::dummy_precision() << std::endl; + std::cout << "Highest = " << Eigen::NumTraits::highest() << std::endl; + std::cout << "Lowest = " << Eigen::NumTraits::lowest() << std::endl; + std::cout << "Digits10 = " << Eigen::NumTraits::digits10() << std::endl; + + example1(); + example2(); + example3(); + example4(); + example5(); + example7(); + example8(); + example9(); + example11(); + example12(); + example13(); + example14(); + example15(); + example16(); + example17(); + example18a(); + example19(); + example20(); + example21(); + example22(); + // example23(); //requires comparisons. + example24(); +} -- cgit v1.2.3