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