Reproducer for Futility eigenvalue solver test failure

[mattbement]

There is a Futility test that seeks to solve a generalized eigenvalue problem: Ax = B \lambda x, with A = [[0.0015, 0.325],[0., 0.]], B=[[0.1208, 0],[-0.117, 0.184]]. Because of a default subspace dimension of 25, and the requirement of Anasazi (at least in 12.12.1) to not have the subspace dimension be larger than the problem dimension, A and B are expanded to a dimension of 30, with A just adding zeros, and B adding 0.1 on the diagonal. This test fails in setup_solver with the exception below. The reproducer is attached

*** ForTrilinos caught exception!
In ForTrilinos::TrilinosEigenSolver::setup_solver(Teuchos::RCP< Teuchos::ParameterList > const &): /usr/projects/shavanodev/bement/ft/Futility/Trilinos/packages/teuchos/core/src/Teuchos_ArrayView.hpp:568:

Throw number = 28

Throw test that evaluated to true: !( ( 0 <= offset && offset+size_in <= this->size() ) && size_in >= 0 )

Teuchos::ArrayView<long long const>::assert_in_range(): Error, [offset,offset+size) = [2,3) does not lie in the range [0,2)!
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG
STOP 1

futility_ev_failure_reproducer.f90.txt

[aprokop]

@mattbement Thanks a lot! I'll work on this.

[aprokop]

@mattbement There are at least two issues here.

1. Column map Because the values of A for row 2 (and rows 3-30, optional) are not set, the constructed column map by default only has two indices in it coming from the first row column indices. While this is perfectly fine when trying to do matrix-vector multiplications, unfortunately many preconditioners in Ifpack2 expect row map to be a subset of the column map. Moreover, they expect them to be compatible, i.e. row map indices are in the beginning of the column map and in the same order. For serial runs, the fix here is trivial through using a slightly different TpetraCrsMatrix constructor that takes in two maps: row and (new) column. The column map can be the same a row map here.
2. Once I changed 1, and set the diagonal of A to be 0, this caused a different issue. The proposed solver parameters use MueLu with RILUK as a smoother. This is problematic as RILUK expects nonzero diagonal for factorization. Thus, I had to change the diagonal values of A to a nonzero (I set it to 1e-14). This passed the errors for preconditioners.
3. Once 2 was fixed, I started encountering floating point exceptions stemming from LAPACK as called by Anasazi. I'm not sure what the problem is, yet.

[aprokop]

Opened an issue to keep track of possibility of alternative approach: trilinos/trilinos#2544.

[aprokop]

I think it works now. And I currently have the encompassing feeling of rage... I spent several hours trying to figure out the LAPACK FPE exception only to find that it completely disappeared after rebuilding from scratch... :hair_pulling:

The good news is that it seems to work even without setting additional diagonal elements of A and B. I changed the preconditioner to IDENTITY, though, to allow that. If one works with larger matrices and needs relaxation, this may require setting it back.

[aprokop]

Closing for now as #227 was merged. Reopen if necessary.