1 #include "main.h" 2 #include <Eigen/MPRealSupport> 3 #include <Eigen/LU> 4 #include <Eigen/Eigenvalues> 5 #include <sstream> 6 7 using namespace mpfr; 8 using namespace std; 9 using namespace Eigen; 10 11 void test_mpreal_support() 12 { 13 // set precision to 256 bits (double has only 53 bits) 14 mpreal::set_default_prec(256); 15 typedef Matrix<mpreal,Eigen::Dynamic,Eigen::Dynamic> MatrixXmp; 16 17 std::cerr << "epsilon = " << NumTraits<mpreal>::epsilon() << "\n"; 18 std::cerr << "dummy_precision = " << NumTraits<mpreal>::dummy_precision() << "\n"; 19 std::cerr << "highest = " << NumTraits<mpreal>::highest() << "\n"; 20 std::cerr << "lowest = " << NumTraits<mpreal>::lowest() << "\n"; 21 22 for(int i = 0; i < g_repeat; i++) { 23 int s = Eigen::internal::random<int>(1,100); 24 MatrixXmp A = MatrixXmp::Random(s,s); 25 MatrixXmp B = MatrixXmp::Random(s,s); 26 MatrixXmp S = A.adjoint() * A; 27 MatrixXmp X; 28 29 // Basic stuffs 30 VERIFY_IS_APPROX(A.real(), A); 31 VERIFY(Eigen::internal::isApprox(A.array().abs2().sum(), A.squaredNorm())); 32 VERIFY_IS_APPROX(A.array().exp(), exp(A.array())); 33 VERIFY_IS_APPROX(A.array().abs2().sqrt(), A.array().abs()); 34 VERIFY_IS_APPROX(A.array().sin(), sin(A.array())); 35 VERIFY_IS_APPROX(A.array().cos(), cos(A.array())); 36 37 38 // Cholesky 39 X = S.selfadjointView<Lower>().llt().solve(B); 40 VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B); 41 42 // partial LU 43 X = A.lu().solve(B); 44 VERIFY_IS_APPROX((A*X).eval(),B); 45 46 // symmetric eigenvalues 47 SelfAdjointEigenSolver<MatrixXmp> eig(S); 48 VERIFY_IS_EQUAL(eig.info(), Success); 49 VERIFY_IS_APPROX((S.selfadjointView<Lower>() * eig.eigenvectors()), 50 eig.eigenvectors() * eig.eigenvalues().asDiagonal()); 51 } 52 53 { 54 MatrixXmp A(8,3); A.setRandom(); 55 // test output (interesting things happen in this code) 56 std::stringstream stream; 57 stream << A; 58 } 59 } 60 61 extern "C" { 62 #include "mpreal/dlmalloc.c" 63 } 64 #include "mpreal/mpreal.cpp" 65
