GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-13 16:41:50 +08:00
Code Issues Packages Projects Releases Wiki Activity

Extend mpreal unit test to check LLT with complexes.

Browse Source
This commit is contained in:
Gael Guennebaud 2016-11-05 11:28:53 +01:00
parent a354c3ca59
commit 55b4fd1d40
1 changed files with 8 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
unsupported/test/mpreal_support.cpp
View File

@ -12,6 +12,7 @@ void test_mpreal_support()
// set precision to 256 bits (double has only 53 bits) // set precision to 256 bits (double has only 53 bits)
mpreal::set_default_prec(256); mpreal::set_default_prec(256);
typedef Matrix<mpreal,Eigen::Dynamic,Eigen::Dynamic> MatrixXmp; typedef Matrix<mpreal,Eigen::Dynamic,Eigen::Dynamic> MatrixXmp;
typedef Matrix<std::complex<mpreal>,Eigen::Dynamic,Eigen::Dynamic> MatrixXcmp;
std::cerr << "epsilon = " << NumTraits<mpreal>::epsilon() << "\n"; std::cerr << "epsilon = " << NumTraits<mpreal>::epsilon() << "\n";
std::cerr << "dummy_precision = " << NumTraits<mpreal>::dummy_precision() << "\n"; std::cerr << "dummy_precision = " << NumTraits<mpreal>::dummy_precision() << "\n";
@ -25,6 +26,10 @@ void test_mpreal_support()
MatrixXmp B = MatrixXmp::Random(s,s); MatrixXmp B = MatrixXmp::Random(s,s);
MatrixXmp S = A.adjoint() * A; MatrixXmp S = A.adjoint() * A;
MatrixXmp X; MatrixXmp X;
MatrixXcmp Ac = MatrixXcmp::Random(s,s);
MatrixXcmp Bc = MatrixXcmp::Random(s,s);
MatrixXcmp Sc = Ac.adjoint() * Ac;
MatrixXcmp Xc;
// Basic stuffs // Basic stuffs
VERIFY_IS_APPROX(A.real(), A); VERIFY_IS_APPROX(A.real(), A);
@ -37,6 +42,9 @@ void test_mpreal_support()
// Cholesky // Cholesky
X = S.selfadjointView<Lower>().llt().solve(B); X = S.selfadjointView<Lower>().llt().solve(B);
VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B); VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B);
Xc = Sc.selfadjointView<Lower>().llt().solve(Bc);
VERIFY_IS_APPROX((Sc.selfadjointView<Lower>()*Xc).eval(),Bc);
// partial LU // partial LU
X = A.lu().solve(B); X = A.lu().solve(B);