Hi,

currently I am trying to build a solver for finite element analyses (FEA) with the help of Eigen. I am using Eigen version 3.2.6. The matrix which is relevant for solving the problem is very sparse and symmetric, so my first try was to use the ConjugateGradient solver with Upper as specifier (I only store the upper half of the matrix). For small problems, the calculation is very fast and gives the exact same results as other FEA-systems.

Now, for a little larger problem (which is still small in the terms of FEA), i.e. about 60000 x 60000, the solver is extremely slow and also gives completely wrong results. It takes about 8 minutes, the full number of iterations (60000) and has an estimated error of 3.15e+065.

As an alternative I tried the BiCGSTAB solver. For this I had to fill the other half of the matrix. Unfortunately, this didn't change anything. It still takes very long and gives wrong results.

For compiling, I tried different settings. Currently I am using the following flags: EIGEN_VECTORIZE, -O3, -fopenmp. I have no idea if this makes sense or whatsoever. The compiler is GCC (respectively g++) 4.8.1.

Is there something I am doing wrong at compiling? Or is the matrix too big? Is there another (maybe better) way of solving the problem?

Thanks a lot in advance!

Korben


Edit:
Ok, I've added the compiler flag NDEBUG, which made the whole thing a lot faster. Still, ConjugateGradient took about 2 minutes and gave different results than without NDEBUG. BiCGSTAB took 8 minutes and gave wrong results (estimated error: 3.15e+065). I understand that, with NDEBUG defined, things work faster, but how is it possible that results with the same solver differ?

As for the actual FEA: I am now using a direct solver (SparseLU). It takes 2 seconds and still gives (different) wrong results. But this probably just means that I have a mistake somewhere else in the code.

Thanks again.

The problem is that the solver is not converging, so indeed, there might be a problem in the way you assembled the problem. BTW, if your domain is 2D, then the fastest method is SimplicialLDLT<>, if it's 3D, then CG is indeed the best choice.