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