mpreal_support.cpp
Go to the documentation of this file.
00001 #include "main.h"
00002 #include <Eigen/MPRealSupport>
00003 #include <Eigen/LU>
00004 #include <Eigen/Eigenvalues>
00005 
00006 using namespace mpfr;
00007 using namespace std;
00008 using namespace Eigen;
00009 
00010 void test_mpreal_support()
00011 {
00012   // set precision to 256 bits (double has only 53 bits)
00013   mpreal::set_default_prec(256);
00014   typedef Matrix<mpreal,Eigen::Dynamic,Eigen::Dynamic> MatrixXmp;
00015 
00016   std::cerr << "epsilon =         " << NumTraits<mpreal>::epsilon() << "\n";
00017   std::cerr << "dummy_precision = " << NumTraits<mpreal>::dummy_precision() << "\n";
00018   std::cerr << "highest =         " << NumTraits<mpreal>::highest() << "\n";
00019   std::cerr << "lowest =          " << NumTraits<mpreal>::lowest() << "\n";
00020 
00021   for(int i = 0; i < g_repeat; i++) {
00022     int s = Eigen::internal::random<int>(1,100);
00023     MatrixXmp A = MatrixXmp::Random(s,s);
00024     MatrixXmp B = MatrixXmp::Random(s,s);
00025     MatrixXmp S = A.adjoint() * A;
00026     MatrixXmp X;
00027 
00028     // Cholesky
00029     X = S.selfadjointView<Lower>().llt().solve(B);
00030     VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B);
00031 
00032     // partial LU
00033     X = A.lu().solve(B);
00034     VERIFY_IS_APPROX((A*X).eval(),B);
00035 
00036     // symmetric eigenvalues
00037     SelfAdjointEigenSolver<MatrixXmp> eig(S);
00038     VERIFY_IS_EQUAL(eig.info(), Success);
00039     VERIFY_IS_APPROX((S.selfadjointView<Lower>() * eig.eigenvectors()),
00040                       eig.eigenvectors() * eig.eigenvalues().asDiagonal());
00041   }
00042 }
00043 
00044 extern "C" {
00045 #include "mpreal/dlmalloc.c"
00046 }
00047 #include "mpreal/mpreal.cpp"


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:07