Consider Eigen/Suitesparse as a Laspack alternative for "serial" configs #44

Closed
benkirk opened this Issue Feb 17, 2013 · 19 comments

Comments

Projects
None yet
5 participants
Owner

benkirk commented Feb 17, 2013

Both are actively developed and increasingly support threads, which are a significant benefit for shared-memory machines when you don't want to deal with MPI.

Eigen is released under the MPL, so as we discussed before I think we could depend heavily on it and ultimately expect it to be available:
http://eigen.tuxfamily.org/index.php?title=Licensing_FAQ

Suitesparse is GPL, so it shouldn't be a core dependency.

@libmesh-devel, thoughts?

Owner

jwpeterson commented Feb 17, 2013

John

On Feb 17, 2013, at 10:36 AM, "Benjamin S. Kirk" notifications@github.com wrote:

Both are actively developed and increasingly support threads, which are a significant benefit for shared-memory machines when you don't want to deal with MPI.

Eigen is released under the MPL, so as we discussed before I think we could depend heavily on it and ultimately expect it to be available:
http://eigen.tuxfamily.org/index.php?title=Licensing_FAQ

Suitesparse is GPL, so it shouldn't be a core dependency.

@libmesh-devel, thoughts?

Eigen is all headers so its really nice. Also we already configure for it, so it should be easy to just include with libmesh...

On Feb 17, 2013, at 2:52 PM, jwpeterson notifications@github.com wrote:

Eigen is all headers so its really nice. Also we already configure for it, so it should be easy to just include with libmesh…

I've got a branch going for this and will have a pull request soon - you are right, it's easy…

-Ben

Owner

roystgnr commented Feb 17, 2013

Fantastic - thanks, Ben!

Owner

pbauman commented Feb 18, 2013

Rock on, sir. Thanks!

Owner

jwpeterson commented Feb 18, 2013

@benkirk is your plan here to add a new SolverPackage enumeration and LinearSolver subclass? Looks like Eigen supports several sparse iterative and direct solvers...

http://eigen.tuxfamily.org/dox/TutorialSparse.html#TutorialSparseDirectSolvers

Owner

benkirk commented Feb 18, 2013

Precisely. They support some sparse iterative solvers internally (GC and BiCGSTAB are blessed, GMRES is in unsupported but looks like it is incoming), and can interface with SuiteSparse to do sparse direct.

I think we'll easily be able to support the sparse iterative stuff easily enough, as it is all Eigen3, and I'm looking at the sparse direct. I think an option to use an installed SuiteSparse it the way to go (building it is nontrivial), but there is an autotools port that we could fall back to - I'm playing with that now to see how big it would be (https://github.com/benkirk/suitesparse.git)

Owner

pbauman commented Feb 18, 2013

On Mon, Feb 18, 2013 at 10:34 AM, Benjamin S. Kirk <notifications@github.com

wrote:

I think an option to use an installed SuiteSparse it the way to go
(building it is nontrivial),

Check out our build scripts on the ICES system
(/opt/apps/ossw/build_scripts/octave ... I think). I setup a build script
for SuiteSparse that might be helpful.

Owner

jwpeterson commented Feb 18, 2013

Testing gmail filtering... disregard.

Owner

benkirk commented Feb 18, 2013

Make that https://github.com/libMesh/suitesparse.git - I needed to make some changes, forked it to my personal account, but thought it more legitimate to put it somewhere the whole team could commit to, if needed.

This is a pretty slick wrapper:

$ ./build.bash ALL
Owner

benkirk commented Feb 20, 2013

I've got a branch at https://github.com/benkirk/libmesh.git working this, feel free to throw stones...

Owner

benkirk commented Feb 21, 2013

Argh!!! MatType is a template parameter for Eigen, and a #define by PETSc...

let the hackery ensue.

Owner

roystgnr commented Feb 21, 2013

On Thu, 21 Feb 2013, Benjamin S. Kirk wrote:

Argh!!! MatType is a template parameter for Eigen, and a #define by PETSc...

A #define in which PETSc header? It seems like we should be trying to
keep that out of non-petsc-related source files regardless.

Owner

benkirk commented Feb 21, 2013

On Feb 21, 2013, at 10:25 AM, roystgnr notifications@github.com wrote:

A #define in which PETSc header? It seems like we should be trying to
keep that out of non-petsc-related source files regardless.

benkirk(3)$ grep -r MatType /software/x86_64/petsc/3.3-p4/include | grep define
/software/x86_64/petsc/3.3-p4/include/finclude/petscmatdef.h:#define MatType character*(80)

On Thu, 21 Feb 2013, Benjamin S. Kirk wrote:

benkirk(3)$ grep -r MatType /software/x86_64/petsc/3.3-p4/include | grep define
/software/x86_64/petsc/3.3-p4/include/finclude/petscmatdef.h:#define MatType character*(80)

It's also in petscmat.h in their C includes, which gets it into our
petsc_matrix.h, all reasonably enough...

Except that then petsc_matrix.h needs to be included in
sparse_matrix.C for use in SparseMatrix::build(), and so do the
Eigen headers, and then stuff breaks.

I can think of a few ways to fix this, but all are pretty ugly...

Owner

benkirk commented Feb 21, 2013

On Feb 21, 2013, at 10:47 AM, LibMesh Developers notifications@github.com wrote:

I can think of a few ways to fix this, but all are pretty ugly…

For now this seems to do the trick, but it is a one-off…

// hack to avoid MatType collision...
#undef libMeshSaveMatType
#ifdef MatType

define MatType libMeshSaveMatType

undef MatType

#endif

// Eigen includes
#include <Eigen/Dense>
#include <Eigen/Sparse>

#ifdef libMeshSaveMatType

define libMeshSaveMatType MatType

undef libMeshSaveMatType

#endif

On Thu, 21 Feb 2013, Benjamin S. Kirk wrote:

For now this seems to do the trick, but it is a one-off…

Yeah; that'll break just as soon as PETSc replaces that macro with a
typedef, but good enough for now.

I think the long-term solution is going to have to be
{petsc,eigen...}_matrix_build.{h,C} files, each with a
PetscMatrixBuild() type method. Then sparse_matrix.h will only have
to include petsc_matrix_build.h, which won't have to include
petsc_matrix.h, and we'll never have a single source file which
requires two different packages' headers.

benkirk added a commit that referenced this issue Feb 21, 2013

Owner

benkirk commented Feb 21, 2013

Gotta love this one, from Ubuntu-LTS's bundled eigen3:
The sparse module API is not stable yet. To use it anyway, please define the EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET preprocessor token.

  CXX      src/numerics/libmesh_opt_la-eigen_preconditioner.lo
In file included from ./include/libmesh/eigen_core_support.h:45:0,
                 from ./include/libmesh/eigen_preconditioner.h:28,
                 from src/numerics/eigen_preconditioner.C:25:
/usr/include/eigen3/Eigen/Sparse:19:2: error: #error The sparse module API is not stable yet. To use it anyway, please define the EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET preprocessor token.
In file included from src/numerics/eigen_preconditioner.C:26:0:
./include/libmesh/eigen_sparse_matrix.h:89:11: error: 'Triplet' in namespace 'Eigen' does not name a type
make[1]: *** [src/numerics/libmesh_opt_la-eigen_preconditioner.lo] Error 1
make[1]: *** Waiting for unfinished jobs....
make: *** [all-recursive] Error 1
program finished with exit code 2
elapsedTime=483.055488

benkirk added a commit that referenced this issue Feb 21, 2013

Owner

roystgnr commented Apr 17, 2013

Are we ready to close this issue? Eigen isn't behaving quite as well as PETSc yet (see commit f220599 ) but it seems to be failing no more often than Laspack or Trilinos, and presumably for the same PETSc-has-better-default-solvers reasons.

Owner

benkirk commented Apr 17, 2013

No problem - its kinda a reminder to me to update the support as they add incrementally more solvers, but seeing Rhys send out a new announcement whenever they roll a new release accomplishes the same thing!

@roystgnr roystgnr closed this Apr 17, 2013

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment