[136] | 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 |
|
---|
| 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 | }
|
---|