Feature eigen

Context

This MR replaces gmm by Eigen.
Eigen is newer, easier to use, has a larger community and could possibly be faster than gmm. Furthermore, Eigen has an interface to the Intel's Pardiso solver, which might be more efficient than MUMPS.

Tasks

For now, the code has merely been translated, but not really refactored nor optimized.
To-do list:

• replace gmm::matrix<double> by Eigen::Matrix<double>
• replace std::vector<double> by Eigen::(Row)Vector<double>, where applicable and useful
• refactor elementary matrices computation (mainly tbox::Element:build() methods), especially clean the thermo-mechanical matrices
• change method return type, ie void build(Matrix const &A, Matrix &K) should become Matrix K build(Matrix const &A), where applicable (mainly tbox::Element)
• replace tbox::Pt by Eigen::VectorXd
• investigate Eigen multithreading feature
• investigate Pardiso interface

The unchecked tasks will be dealt with in the next MR.

Issues

Eigen uses the method of co-factors to inverse a matrix, whose size is fixed and up to 4x4.For now, the matrices in Mem are dynamically sized, but there should be a way to force Eigen to use the method of co-factors for MemTri3, MemQuad4 and MemTetr4. This could partly address #27.

Notes

I replaced std::vector<double> by Eigen::VectorXd where those vectors were used for linear algebra. There are some exceptions, check the comments in my commits.
The main exception is that some vectors are public and hold variables and results. I decided to keep std::vector for such cases in order to have more flexibility when interfacing in python (SWIG already already fully interfaces std::vector). For linear algebra on the C++ side, these std::vector are mapped with Eigen::Map<Eigen::Vector>, which holds a reference to std::vector to keep it up to date, and can be manipulated as Eigen::Vector.

Additional changes

• Git LFS support
• Force to look for MKL before looking for other BLAS
• Minor changes in flow solver

Tests

Passed on ubuntu18.04 (python 2.7.15), debian4.9 (python 2.7.13), msys2 (python 3.8.1).
@R.Boman and @Kim.Liegeois: Please double check that the solvers you implemented still give accurate results, particularly using multithreading. The battery is passing, but I did not run and check each test individually (only some of them)!

Here we use a Map instead of a Vector, because we want to keep the public variables as std::vector (easier in pyhton)

Here we keep std::vector and not Map, because b is used in Assign::apply(std::vector &), so it is easier.

When using Eigen, I was not able to decrease the residual below 3e-11. 1e-10 is reached at iteration 70 instead of 67, as previously with gmm.

no problem; the solution is the same

This file needs to be double checked and cleaned (@Kim.Liegeois )

This fails and gives weird errors involving... tbb?! Needs to be investigated...

lapacke is not easily available on windows with msvc

Does not explain the error, but I can remove it if it is not readily available on one OS.

Map instead of Vector to be consistent with Solver classes in modules. Is there a template class that superseeds both Map and Vector?

If you want the method to accept both std::vector and Eigen::VectorXd types, couldn't you use Eigen::VectorXd const & as argument?

• The call with a Eigen::VectorXd would be setSoE(A, b)
• The call with a std::vector b would be setSoE(A, Eigen::Map<Eigen::VectorXd>(b)) In this way, the conversion is explicit during the call and it is made outside the implementation of the called class.

If it is possible to do so, it is cleaner because we don't need to know what happens somewhere else in the code when we try to understand the types of the arguments of this method. The method also become compatible with a classical Eigen::VectorXd object (which is the type that we expect here).

What I would have liked is a b as Eigen::Vector or as Eigen::MapEigen::Vector, because b can be:

• the RHS of a Picard iteration soe, which we are not really interested in (thus could be a temporary Vector)
• the RHS of a Newton iteration soe, ie, the residual, which were are interested in (thus has to be a Map)

However, I see no way of doing that...

The question will be solved if we find an easy way to wrap Eigen::Vectors in python.

I'll look into it... in the future!

Switch to array for coefficient-wise division.

Map instead of Vector to be consistent with Solver classes in modules. Is there a template class that superseeds both Map and Vector?

It looks that it is rather simple to convert a std::vector into an Eigen::VectorXd using a Map. The opposite is impossible (std::vector cannot "share" its memory). Thus, I don't see any (clean) way to accept both std::vector and Eigen::VectorXd in this case.

At this point, I have noticed no difference in performance between Eigen and gmm.

added cleaning enhancement labels

I could not add reviewers (maybe because of the WIP status).

marked the checklist item refactor elementary matrices computation (mainly tbox::Element:build() methods), especially clean the thermo-mechanical matrices as completed

marked the checklist item refactor elementary matrices computation (mainly tbox::Element:build() methods), especially clean the thermo-mechanical matrices as incomplete

added 1 commit

• 83e47124 - fix for tbb

Compare with previous version

I am not sure why I need to block here... All idxUp and idxLw are different, so there should not be any thread conflicts... @R.Boman any ideas what I am missing?

I guess that the matrix is structurally modified by this code. How does adding a zero is handled internally? I don't know, and even if it does not modify the internal structure of the matrix today, can you make this assumption for the future version of Eigen?

Check the structure of the matrix before and after this loop. It should be the same.

You could also add a debug loop before that, which would check whether all the indices that you access have been previously allocated (and are still allocated after the loop).

Note: I would add : assert(tUp.size()<= pbl->msh->nodes.size() / 10) after the loop to make the assumption explicit. Because, as you know, this could be also a problem.

I guess that the matrix is structurally modified by this code.

I came up with the same conclusion, but only existing coefficients are set to 0, which should not modify the matrix structure. Only the coeffRef method might insert elements and modify the memory.

Note: I would add : assert(tUp.size()<= pbl->msh->nodes.size() / 10) after the loop to make the assumption explicit. Because, as you know, this could be also a problem.

I can add that. However, I doubt that this is the issue since the number of entries in tUp will be equal to the the number of nodes that are connected to the node at ìdxUp` which will be much smaller than the problem size.

Should be solver_ref here, but that does not explain why the pipeline failed on this test...

meshDef has a long history of repeated failures

Okay, that's gmsh related. It needs 8 iterations instead of 5 with gmsh3 -_-

changed the description

added 3 commits

• 3c14d4b0 - Added iteration counter to tests, to increase robustness
• 10a8050f - Corrected mistake
• 0cf68e6c - Merge branch 'adrien' into feature_eigen

Compare with previous version

added 1 commit

• ebaa8863 - Removed parallel loop and increased test tolerance for gmsh3.

Compare with previous version

changed the description