Hey everyone,

i'm kinda new to eigen and i'm currently exploring some of it's features. I'd like to solve a standard task A*x = b with A being sparse and relatively high dimensional. Yet, i've implemented a gauß-Jordan algorithm and now i'm trying compare eigen's sparse lu solver to it. So here is the important part of my code:

Code: Select all

   typedef Eigen::Triplet<double> T;
   std::vector<T> tripletList;
   tripletList.reserve(mxBig.GetColumnCount() * 5);
   for (int i = 0; i < mxBig.GetColumnCount(); i++) {
      for (int j = 0; j < mxBig.GetRowCount(); i++) {
         if (!Null(mxBig.at(i, j)))
            tripletList.push_back(T(i, j, mxBig.at(i, j)));
      }
   }
   Eigen::SparseMatrix<double> A(mxBig.GetColumnCount(), mxBig.GetRowCount());
   A.setFromTriplets(tripletList.begin(), tripletList.end());
   
   Eigen::VectorXd b(vctBigVM.size());
   for (int i = 0; i < vctBigVM.size(); i++) {
      b(i) = vctBigVM.at(i);
   }
   
   Eigen::VectorXd x(vctBigVM.size());
   Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
   solver.compute(A);
   if (solver.info() != Eigen::Success) {
      return false;
   }
   x = solver.solve(b);
   if (solver.info() != Eigen::Success) {
      return false;
   }


mxBig, resp. vctBigVM represent the matrix A and the vector b. It seems to me that the code fragment is relatively close to the tutorial on the eigen website, though i'm getting the following error:

Code: Select all

error C2220: warning treated as error, no object file generated.
warining 2 warning C4714: "const Eigen::Matrix<double,-1,-1,0,-1,-1> Eigen::DenseBase<Derived>::eval(void) const", marked as __forceinline, not inlined ...


which is caused by the line

Code: Select all

x = solver.solve(b);

Can anyone tell me what i'm doing wrong? I'm totally clueless at the moment..

Thanks in advance,
Frederic

This is a warning telling that it choses not to inline a function we asked him to inline. Just remove the flag to treat warnings as warnings, and not as error.

Thanks for this hint. My bad..

By now I've got some results, but they are suprising me quite a bit. My example uses a 50x50 banded matrix having a maximum bandwidth of 5. The solution i get from the self-implemented Gaus-Jordan-Algorithm looks like

Code: Select all

0.00000000000000000 0.00044349400879672946 0.011969600962114784 0.00051326991605071435 0.022737250242648110 0.00058524608827432685 0.031649714464784318 0.00059219085651881207 0.040079377222404136 0.00051441588160154085 0.046464891854909096 0.00031318109509990245 0.048663810675365356 -5,02535714E+11 0.040712896996130517 -0.00082891550100142194 0.00000000000000000 -0.0025011575352028707 -0.19436114335514806 -0.0044055888137627969 -0.35311393992250817 -0.0034694500371233498 -0.43484049276146403 -0.0018905294749781614 -0.46303868908470114 3,90930208E+11 -0.43258873476241727 0.0019593331916156944 -0.34917336342417660 0.0035101067787728195 -0.18999245136724766 0.0043816013361896123 0.00000000000000000 0.0023447854519356376 0.036690150522439145 0.00060310206842460149 0.040626251258077374 -0.00013805353602555447 0.035910452818890992 -0.00043595031769130973 0.028052568874538556 -0.00055176113841027804 0.019490041339801065 -0.00053514928591078144 0.011890531169670098 -0.00043577804792113928 0.0048324997015697178 -0.00026454143389605518 0.00000000000000000 -0.00013774089829451415

,whereas the solution i get from f.e. eigen's lu solver looks like

Code: Select all

0.00000000000000000 5,09275473E+77 1,26112140E+79 4,65121289E+77 2,08548306E+79 3,71771095E+77 2,56674838E+79 2,64823246E+77 2,86467453E+79 1,27318868E+77 2,93342672E+79 -4,07420378E+76 2,72717016E+79 -2,39359472E+77 1,97252322E+79 -5,43244147E+77 0.00000000000000000 -1,01855095E+78 -7,24336394E+79 -1,70423944E+78 -1,41171161E+80 -1,75190763E+78 -1,90672737E+80 -1,50745540E+78 -2,29173963E+80 -1,01855095E+78 -2,49341272E+80 -2,85194265E+77 -2,43841096E+80 6,92614643E+77 -1,86257930E+80 2,32922230E+78 0.00000000000000000 5,09275473E+78 1,47925229E+80 5,85090803E+78 2,61749514E+80 5,28067738E+78 3,35706243E+80 4,16383627E+78 3,87781443E+80 2,45725416E+78 4,08834127E+80 1,60931049E+77 3,89723310E+80 -2,72513306E+78 2,89802464E+80 -7,35867918E+78 0.00000000000000000 -1,49217714E+79

I've tried different solvers as well as the rather naive call of the inverse() method, but i'm always getting those extremly high values. Can somebody provide a guess or maybe even give an explanation on whats going wrong?

Thanks in advance,
Frederic

Hi Frederic,

I am totally new to Eigen but maybe there is an error in the code you posted:

Code: Select all

   typedef Eigen::Triplet<double> T;
   std::vector<T> tripletList;
   tripletList.reserve(mxBig.GetColumnCount() * 5);
   for (int i = 0; i < mxBig.GetColumnCount(); i++) {
      for (int j = 0; j < mxBig.GetRowCount(); i++) {
         if (!Null(mxBig.at(i, j)))
            tripletList.push_back(T(i, j, mxBig.at(i, j)));
   }


It seems to me that in the second loop you are using "i++" instead of "j++"; maybe this is the cause of the weird results..

Pietro

Hey Pietro,

thanks for your reply, but unfortunately that was just a typo in the board entry. Of course it needs to be j++ and my code already says so.. Further more I've checked the input of the solvers more than twice and i'm relatively sure they are correct.

Regards,
Frederic

Can you post your matrix in the matrix-market format (don't forget to zip it!):

Code: Select all

#include <unsupported/Eigen/SparseExtra>
...
Eigen::saveMarket(spmat, "filename.mtx");


Hey!

Of course I can. Here it is: A & b as .mtx (zipped)

From the data you provided I get:

Code: Select all

0  5.09275e+61  1.26112e+63  4.65121e+61  2.08548e+63  3.71771e+61  2.56675e+63  2.64823e+61  2.86467e+63  1.27319e+61  2.93343e+63  -4.0742e+60  2.72717e+63 -2.39359e+61  1.97252e+63 -5.43244e+61            0 -1.01855e+62 -7.24336e+63 -1.70424e+62 -1.41171e+64 -1.75191e+62 -1.90673e+64 -1.50746e+62 -2.29174e+64 -1.01855e+62 -2.49341e+64 -2.85194e+61 -2.43841e+64  6.92615e+61 -1.86258e+64  2.32922e+62            0  5.09275e+62  1.47925e+64  5.85091e+62   2.6175e+64  5.28068e+62  3.35706e+64  4.16384e+62  3.87781e+64  2.45725e+62  4.08834e+64  1.60931e+61  3.89723e+64 -2.72513e+62  2.89802e+64 -7.35868e+62            0 -1.49218e+63


with all solvers including non-Eigen ones with a relative error in the order of 1e-16. Makes sure that your Gauss-Jordan implementation compute a correct solution by checking |A*x-b|/|b|. Then makes sure you are feeding both solver with the exact same values.

BTW, you might also be interested by our benchmark routine: http://eigen.tuxfamily.org/dox/group__T ... tml#title3

Thanks for your input ggael,

i finally got it. There was a single value being copied wrong, which caused all this drama. I've one last question regarding the overview of the solvers. It says that SimplicialLDLT is recommended for "not too large problems". Is there a more precise answer on what "not too large" means? dim = 10, 100, 1.000, 10.000, ...?

"not too large problems" should rather be read as not too dense Typically, it will work well for FEM on a 2D domain, but lack performance for 3D domains.