00001 #include "main.h" 00002 #include <Eigen/MPRealSupport> 00003 #include <Eigen/LU> 00004 #include <Eigen/Eigenvalues> 00005 00006 using namespace mpfr; 00007 using namespace std; 00008 using namespace Eigen; 00009 00010 void test_mpreal_support() 00011 { 00012 // set precision to 256 bits (double has only 53 bits) 00013 mpreal::set_default_prec(256); 00014 typedef Matrix<mpreal,Eigen::Dynamic,Eigen::Dynamic> MatrixXmp; 00015 00016 std::cerr << "epsilon = " << NumTraits<mpreal>::epsilon() << "\n"; 00017 std::cerr << "dummy_precision = " << NumTraits<mpreal>::dummy_precision() << "\n"; 00018 std::cerr << "highest = " << NumTraits<mpreal>::highest() << "\n"; 00019 std::cerr << "lowest = " << NumTraits<mpreal>::lowest() << "\n"; 00020 00021 for(int i = 0; i < g_repeat; i++) { 00022 int s = Eigen::internal::random<int>(1,100); 00023 MatrixXmp A = MatrixXmp::Random(s,s); 00024 MatrixXmp B = MatrixXmp::Random(s,s); 00025 MatrixXmp S = A.adjoint() * A; 00026 MatrixXmp X; 00027 00028 // Cholesky 00029 X = S.selfadjointView<Lower>().llt().solve(B); 00030 VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B); 00031 00032 // partial LU 00033 X = A.lu().solve(B); 00034 VERIFY_IS_APPROX((A*X).eval(),B); 00035 00036 // symmetric eigenvalues 00037 SelfAdjointEigenSolver<MatrixXmp> eig(S); 00038 VERIFY_IS_EQUAL(eig.info(), Success); 00039 VERIFY_IS_APPROX((S.selfadjointView<Lower>() * eig.eigenvectors()), 00040 eig.eigenvectors() * eig.eigenvalues().asDiagonal()); 00041 } 00042 } 00043 00044 #include "mpreal.cpp"