| [eigen] boost::multiprecision : wrong eigen values |
[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]
Dear All, in case you are interested, I've attached an example where the eigen values using boost::multiprecision are wrong. In addition, in contrast to the documentation, they are not sorted, which could be a hint for the problem. I included a Makefile, which provides executables for double, long double, boost::multiprecision and mpfr. All, except boost::multiprecision, are consistent. As you can also realize, boost::multiprecision is much slower compared to mpfr, so I'm not using it and I have no need for it to be fixed. However, it might be interesting to see, why it fails. In case you wonder how I come up with the matrix, it's just an example I found in an application of mine. Best regards, Peter
#include <stdlib.h>
#include <iostream>
#include <cmath>
#include <Eigen/Dense>
#include <algorithm>
// --------------------
#if defined( D_boost )
// --------------------
#include <boost/multiprecision/cpp_bin_float.hpp>
#define MyBoost_Prec 500
const int MyBoostPrec2 = MyBoost_Prec * log(2.0)/log(10.0) ;
typedef boost::multiprecision::cpp_bin_float< MyBoostPrec2 > mp_backend; ///< Why is it using decimal precision???
typedef boost::multiprecision::number<mp_backend, boost::multiprecision::et_off> TpFloat;
// courtesy of Carsten Daniel Burgard
// we have to turn off expression templates to make Eigen like boost
// this is apparently needed because multiprecision doesn't come with an int() operator
namespace Eigen {
namespace internal {
template<typename NewType>
struct cast_impl<TpFloat,NewType> {
static inline NewType run(const TpFloat& x) {
return x.convert_to<NewType>();
}
};
}
}
// --------------------
#elif defined( D_mpfr )
// --------------------
#include <unsupported/Eigen/MPRealSupport>
typedef mpfr::mpreal TpFloat;
// --------------------
#elif defined( D_ld )
// --------------------
typedef long double TpFloat;
// --------------------
#elif defined ( D_d )
// --------------------
typedef double TpFloat;
#endif
// ------------------------------------------------------------
typedef Eigen::Matrix< TpFloat, Eigen::Dynamic, 1> TpVector;
typedef Eigen::Matrix< TpFloat, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> TpMatrix;
typedef Eigen::Array< TpFloat , Eigen::Dynamic, 1> TpArray;
int main( const int argc, const char *argv[] )
{
const int dim = 32;
TpMatrix H( dim, dim );
TpArray E,EE;
int Precision = (argc>1) ? atoi(argv[1]) : 500 ;
#if defined( D_boost )
Precision = MyBoost_Prec;
#elif defined( D_mpfr )
mpfr::mpreal::set_default_prec(Precision);
#else
Precision = 8 * sizeof( TpFloat );
#endif
H << -0.125, -0, -0, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0, -0.7071067812, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, -0.125, -0, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0, 0, 0, -0.7071067812, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, -0, -0.125, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0, 0, 0, -0.7071067812, 0, 0, 0, 0, 0, 0, 0,
-0, -0, -0, -0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0, -0.7071067812, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, -0.125, -0, -0, -0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, -0, -0.125, -0, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, -0, -0, -0.125, -0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, -0, -0, -0, -0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, -0.125, -0, -0, -0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, -0, -0.125, -0, -0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, -0, -0, -0.125, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, -0, -0, -0, -0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.7071067812, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0, 0,
0, -0.7071067812, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.7071067812, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7071067812, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, -0.7071067812, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7071067812, 0,
0, -0.7071067812, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0,
0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, -0.7071067812, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, -0.7071067812, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7071067812, 0,
0, 0, 0, -0.7071067812, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7071067812,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0, 0.7071067812, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.375, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0, 0, 0, -0.7071067812, 0, 0, 0, 0, 0, 0, 0, 0, 0.375, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7071067812, 0, 0, 0, 0, 0.7071067812, 0, 0, 0, 0, 0, 0.375, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7071067812, 0, 0, 0.7071067812, 0, 0, 0, 0, 0, 0.375;
Eigen::SelfAdjointEigenSolver<TpMatrix> eigensolverH;
eigensolverH.compute( H );
EE =eigensolverH.eigenvalues().array();
std::sort( EE.data(), EE.data() + dim ); // needed for mp-boost, should be sorted according to docs.
E = EE - EE(0);
std::cout << E.transpose() << std::endl;
return 0;
}
CFLAGS=-O3 -march=native -DNDEBUG All: MpTest_d MpTest_ld MpTest_mpfr MpTest_boost clean: rm -f MpTest_d MpTest_ld MpTest_mpfr MpTest_boost MpTest_d: MpTest.C g++ -DD_d $(CFLAGS) MpTest.C -I $(HOME)/Eigen/eigen -o MpTest_d MpTest_ld: MpTest.C g++ -DD_ld $(CFLAGS) MpTest.C -I $(HOME)/Eigen/eigen -o MpTest_ld MpTest_mpfr: MpTest.C g++ -DD_mpfr $(CFLAGS) MpTest.C -I $(HOME)/Eigen/eigen -I $(HOME)/Eigen/eigen/unsupported/test/mpreal -l mpfr -o MpTest_mpfr MpTest_boost: MpTest.C g++ -DD_boost $(CFLAGS) MpTest.C -I $(HOME)/Eigen/eigen -I $(HOME)/RG/Boost/boost-stable -o MpTest_boost
| Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |