mpreal_support.cpp
Go to the documentation of this file.
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 
11 EIGEN_DECLARE_TEST(mpreal_support)
12 {
13  // set precision to 256 bits (double has only 53 bits)
14  mpreal::set_default_prec(256);
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
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 }
Jet< T, N > cos(const Jet< T, N > &f)
Definition: jet.h:426
EIGEN_DEVICE_FUNC ComputationInfo info() const
Reports whether previous computation was successful.
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
DiscreteKey S(1, 2)
static double epsilon
Definition: testRot3.cpp:37
EIGEN_DEVICE_FUNC const ExpReturnType exp() const
#define VERIFY_IS_APPROX(a, b)
EIGEN_DONT_INLINE void llt(const Mat &A, const Mat &B, Mat &C)
Definition: llt.cpp:5
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
EIGEN_DECLARE_TEST(mpreal_support)
EIGEN_DEVICE_FUNC const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
static int g_repeat
Definition: main.h:169
SelfAdjointEigenSolver< PlainMatrixType > eig(mat, computeVectors?ComputeEigenvectors:EigenvaluesOnly)
RealScalar s
#define VERIFY(a)
Definition: main.h:380
internal::nested_eval< T, 1 >::type eval(const T &xpr)
EIGEN_DEVICE_FUNC const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
const int Dynamic
Definition: Constants.h:22
The matrix class, also used for vectors and row-vectors.
EIGEN_DEVICE_FUNC bool isApprox(const Scalar &x, const Scalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
#define X
Definition: icosphere.cpp:20


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:56