Unexpected segmentation fault with hypre and SuperLU

[andrea-franceschini]

@castelletto1 has pointed my attention to the following issue.
Compiling geosx using hypre linear algebra interface and running the integrated test compositional_multiphase_wells_1d.xml, at the first linear system solution, we have:

stack trace

Received signal 8: Floating point exception
Floating point exception:

** StackTrace of 23 frames **
Frame 1: cxx_utilities::handler(int, int, int)
Frame 2: cxx_utilities::handler1(int)
Frame 3:
Frame 4: dtrsm_ounncopy_HASWELL
Frame 5: dtrsm_RNUN
Frame 6: dtrsm_
Frame 7: pdgstrf2_trsm
Frame 8: pdgstrf
Frame 9: pdgssvx
Frame 10: hypre_SLUDistSetup
Frame 11: geosx::HypreSolver::solve_direct(geosx::HypreMatrix&, geosx::HypreVector&, geosx::HypreVector&)
Frame 12: geosx::HypreSolver::solve(geosx::HypreMatrix&, geosx::HypreVector&, geosx::HypreVector&)
Frame 13: geosx::SolverBase::SolveSystem(geosx::DofManager const&, geosx::HypreMatrix&, geosx::HypreVector&, geosx::HypreVector&)
Frame 14: geosx::ReservoirSolver::SolveSystem(geosx::DofManager const&, geosx::HypreMatrix&, geosx::HypreVector&, geosx::HypreVector&)
Frame 15: geosx::SolverBase::NonlinearImplicitStep(double const&, double const&, int, geosx::DomainPartition*, geosx::DofManager const&, geosx::HypreMatrix&, geosx::HypreVector&, geosx::HypreVector&)
Frame 16: geosx::ReservoirSolver::SolverStep(double const&, double const&, int, geosx::DomainPartition*)
Frame 17: geosx::SolverBase::Execute(double, double, int, int, double, geosx::dataRepository::Group*)
Frame 18: geosx::EventBase::Execute(double, double, int, int, double, geosx::dataRepository::Group*)
Frame 19: geosx::EventManager::Run(geosx::dataRepository::Group*)
Frame 20: geosx::ProblemManager::RunSimulation()
Frame 21: main
Frame 22: __libc_start_main
Frame 23: _start

It means that the BLAS routine dtrsm throughs an error when solving a supernode (to be precise, the first one!). Thus, it seems that the matrix cannot be solved and the only reason for this is that it's singular. So, I saved the matrix (145 x 145) and computed all its eigenvalues. It seems to be a regular matrix, indeed, both MATLAB and geosx itself, when compiled with Trilinos, are able to solve it.

Moreover, using the saved version of the matrix (all entries are printed with 15 digits) and calling directly either SuperLU or hypre (that calls SuperLU), everything is fine. The matrix doesn't show any problem!

Digging into the problem, it seems that the geosx + hypre + SuperLU version starts to be different from the hypre + SuperLU and SuperLU versions already at the permuting stage (before the factorization). That stage, based on both pattern and numerical entries, produces different results and from that moment, the two factorizations are very different (with one producing a wrong supernode, that the BLAS routine dtrsm is not able to solve).

Finally, I tried a very simple trick: to add 1.e-20 to the diagonal terms of the matrix. With this simple change, the geosx + hypre + SuperLU produces the same exact result as the two drivers that don't use geosx.