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
EIGEN_DECLARE_TEST(mpreal_support)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 }
67