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
test_mpreal_support()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