mpreal_support.cpp (2439B)
1 #include <mpreal.h> // Must be included before main.h. 2 #include "main.h" 3 #include <Eigen/MPRealSupport> 4 #include <Eigen/LU> 5 #include <Eigen/Eigenvalues> 6 #include <sstream> 7 8 using namespace mpfr; 9 using namespace Eigen; 10 11 EIGEN_DECLARE_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 typedef Matrix<std::complex<mpreal>,Eigen::Dynamic,Eigen::Dynamic> MatrixXcmp; 17 18 std::cerr << "epsilon = " << NumTraits<mpreal>::epsilon() << "\n"; 19 std::cerr << "dummy_precision = " << NumTraits<mpreal>::dummy_precision() << "\n"; 20 std::cerr << "highest = " << NumTraits<mpreal>::highest() << "\n"; 21 std::cerr << "lowest = " << NumTraits<mpreal>::lowest() << "\n"; 22 std::cerr << "digits10 = " << NumTraits<mpreal>::digits10() << "\n"; 23 24 for(int i = 0; i < g_repeat; i++) { 25 int s = Eigen::internal::random<int>(1,100); 26 MatrixXmp A = MatrixXmp::Random(s,s); 27 MatrixXmp B = MatrixXmp::Random(s,s); 28 MatrixXmp S = A.adjoint() * A; 29 MatrixXmp X; 30 MatrixXcmp Ac = MatrixXcmp::Random(s,s); 31 MatrixXcmp Bc = MatrixXcmp::Random(s,s); 32 MatrixXcmp Sc = Ac.adjoint() * Ac; 33 MatrixXcmp Xc; 34 35 // Basic stuffs 36 VERIFY_IS_APPROX(A.real(), A); 37 VERIFY(Eigen::internal::isApprox(A.array().abs2().sum(), A.squaredNorm())); 38 VERIFY_IS_APPROX(A.array().exp(), exp(A.array())); 39 VERIFY_IS_APPROX(A.array().abs2().sqrt(), A.array().abs()); 40 VERIFY_IS_APPROX(A.array().sin(), sin(A.array())); 41 VERIFY_IS_APPROX(A.array().cos(), cos(A.array())); 42 43 // Cholesky 44 X = S.selfadjointView<Lower>().llt().solve(B); 45 VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B); 46 47 Xc = Sc.selfadjointView<Lower>().llt().solve(Bc); 48 VERIFY_IS_APPROX((Sc.selfadjointView<Lower>()*Xc).eval(),Bc); 49 50 // partial LU 51 X = A.lu().solve(B); 52 VERIFY_IS_APPROX((A*X).eval(),B); 53 54 // symmetric eigenvalues 55 SelfAdjointEigenSolver<MatrixXmp> eig(S); 56 VERIFY_IS_EQUAL(eig.info(), Success); 57 VERIFY( (S.selfadjointView<Lower>() * eig.eigenvectors()).isApprox(eig.eigenvectors() * eig.eigenvalues().asDiagonal(), NumTraits<mpreal>::dummy_precision()*1e3) ); 58 } 59 60 { 61 MatrixXmp A(8,3); A.setRandom(); 62 // test output (interesting things happen in this code) 63 std::stringstream stream; 64 stream << A; 65 } 66 }