Skip to content

SparseModule and MPFR C++.

Submitted by ber..@..ign.fr

Assigned to Nobody

Link to original bugzilla bug (#116)
Version: 3.0
Platform: x86 - 32-bit

Description

Hello
I have tried to use the SparseModule with MPFR C++, but I have an error that I cannot understand. I don't know whether the problem comes from Eigen or from MPFR C++.

My code runs correctly when the numerical type of the matrix is double but does not work when the type is mpreal in Sparse Module.
Can you give me some hint about this error ?

Thank you for your answer.

Bertrand

To reproduce my bug, you can use this code

////////////////////////////////////////////////////////////////////////

#include <Eigen/MPRealSupport>

#include "mpreal.h"

#include <Eigen/Core>

#include <Eigen/LU>

#include <Eigen/Dense>

#include <Eigen/Sparse>

#include <Eigen/SparseExtra>

using namespace mpfr;

using namespace std;

template<class U>

void Test(Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> A ,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> B){

Eigen::Matrix<U,Eigen::Dynamic,Eigen::Dynamic> T_A_;  

Eigen::Matrix<U,Eigen::Dynamic,Eigen::Dynamic> T_B_;  

T_A_.resize(A.rows (),A.cols());   

T_B_.resize(B.rows (),B.cols());   



for(unsigned int r=0;r<A.rows ();r=r+1)  

    for(unsigned int c=0;c<A.cols();c=c+1)  

        T_A_.coeffRef(r,c)=A.coeffRef(r,c);   



for(unsigned int r=0;r<B.rows ();r=r+1)  

    for(unsigned int c=0;c<B.cols();c=c+1)  

        T_B_.coeffRef(r,c)=B.coeffRef(r,c);   





Eigen::Matrix<U,Eigen::Dynamic,Eigen::Dynamic> T_X_;  

  

Eigen::LDLT< Eigen::Matrix<U,Eigen::Dynamic,Eigen::Dynamic> > ldlt(T_A_);  

T_X_=T_B_;  

ldlt.solveInPlace(T_X_);  

std::cout << "relative error: " << (T_A_*T_X_ - T_B_).norm() / T_B_.norm() << std::endl;      

  

  

Eigen::DynamicSparseMatrix&lt;U&gt; T_A ;  

Eigen::Matrix<U,Eigen::Dynamic,Eigen::Dynamic> T_B(B.rows (),B.cols()) ;  

T_A.resize(A.rows (),A.cols());  

for(unsigned int i=0;i<100;i=i+1){  

    for(unsigned int j=0;j<100;j=j+1)  

        T_A.coeffRef(i,j)=T_A_.coeffRef(i,j);   

    T_B.coeffRef(i,0)=T_B_.coeffRef(i,0);  

}

Eigen::Matrix<U,Eigen::Dynamic,Eigen::Dynamic> T_X=T_B;

Eigen::SparseLDLT< Eigen::SparseMatrix<U,Eigen::Upper|Eigen::SelfAdjoint> > llt(T_A);  

if (llt.succeeded()){  

     llt.solveInPlace(T_X);  

  }  

else{  

    std::cout<<"error"<<std::endl;  

    }  

std::cout << "relative error: " << (T_A*T_X - T_B).norm() / T_B.norm() << std::endl;

};

int main(int argc,char**argv)

{

Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> A = Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>::Random(100,100);  

Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> B = Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>::Random(100,1);  

  

mpreal::set_default_prec(256);  

  

Test&lt;double&gt;(A,B);  

Test&lt;mpreal&gt;(A,B);  



return(0);

}
/////////////////////////////////////////////////////////////////////////
relative error: 18.1509
relative error: 22.196
relative error: 18.1509
Error
////////////////////////////////////////////////////////////////////////

Edited by Eigen Bugzilla