Anasazi: mutually exclusive checks

[aprokop]

@trilinos/anasazi

Expectations

I'm trying to compute a single eigenpair of a generalized eigensystem using generalized Davidson method. I'm unable to set the parameters to get it done for 2x2 size matrices.

Current Behavior

Anasazi contains the following block of checks:

TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV()+2 > pl.get<int>("Maximum Subspace Dimension"),
                                std::invalid_argument, "Maximum Subspace Dimension must be strictly greater than NEV");
TEUCHOS_TEST_FOR_EXCEPTION( d_maxSubspaceDim > MVT::GetGlobalLength(*problem->getInitVec()),
                                std::invalid_argument, "Maximum Subspace Dimension cannot exceed problem size");

First of all, in the parameter list I specify "Maximum Subspace Dimension" to be 1 or 2, otherwize the last check fails (remember, the matrices are 2x2). However, setting it this way comes into conflict with the first check where I request a single eigenpair which results in (1+2 > 2).

Motivation and Context

The application in question has to solve a bunch of eigensystems of different sizes, including the small ones. What they typically did was to pad the matrices to 30x30 so that the subspace checks pass, however that leads to a bunch of other hacks (like setting nonzero diagonal for RILUK preconditioners to construct, etc). I would like to directly pass the proper eigensystem to Anasazi. I see this as a reasonable use case.

Possible Solution

I would like somebody to take a look at the checks and double check whether they make sense for smaller sized matrices. Specifically, I'm not sure why +2 is necessary in the check. If there are limitations in Anasazi for those situations, it would be nice to have that documented.

Steps to Reproduce

So far, I have a small standalone Fortran example to reproduce. I don't think it's necessary to put it here, though, as the situation should be clear just by looking at the checks.

Additional Information

The issues is in both develop and 12.12 branch.

[mhoemmen]

@hkthorn

[hkthorn]

@aprokop @mhoemmen Sorry, I was on vacation. I'll look into this. Unfortunately, generalized Davidson was not written by me or Chris Baker, so I am not as familiar with how Steve Hamilton implemented the solver.

[aprokop]

@hkthorn Can you ping Steven on this?

[mhoemmen]

@hkthorn never apologize for vacation :D

[hkthorn]

@aprokop @mhoemmen I'll take a look first to see if I can rationalize the logic. I will pass questions on to Steven for clarification, if necessary.

[hkthorn]

@aprokop @mhoemmen I have taken a look at the Generalized Davidson solver. While I do not know why there is a +2, I can certainly see why there would be a +1. If this is a non-Hermitian generalized eigensystem then there can be complex eigenvalues, which need to be represented by a 2x2 blocks since real arithmetic is being used in this solver. That means there could be an extra Ritz value and vector that need to be stored in the event of a restart.

The check that is failing ensures that the maximum subspace dimension is strictly greater than the number of eigenvalues that are being computed. If that is changed to only include 1 extra eigenvalue to account for potential conjugate pairs then the logic becomes:

d_problem->getNEV()+1 > pl.get("Maximum Subspace Dimension")

This will still throw an error for a 2x2 system if getNEV() returns 1 and the "Maximum Subspace Dimension" is 2. You can't set the "Maximum Subspace Dimension" larger than 2 for obvious reasons (numerically), but also because the next check in the solver would fail:

d_maxSubspaceDim > MVT::GetGlobalLength(*problem->getInitVec())

So, I can discuss with Steven Hamilton about changing inequality to express exactly what it states. I can imagine he didn't expect the solver be used on small systems. Maybe he won't care whether the logic is changed.

[hkthorn]

@aprokop @mhoemmen

Here is Steven's reply to my inquiry:

_Heidi,

You are correct that the issue is due to handling of conjugate pairs, and I might be handling it in an inelegant way.

The problem is that according to the current implementation, the subspace may be expanded by two vectors in a single iteration. If the current “best" Ritz value is a conjugate pair, then we store the corresponding residual as two vectors (effectively the real and imaginary components of the actual complex residual vector). The subspace gets expanded by applying the preconditioner to the residual, which adds two vectors to the current subspace. If you want 10 values/vectors and only allow a subspace size of 11, you’re potentially unable to expand the subspace in the “correct” way.

At the very least, the error message should be corrected since the actual check is more stringent than what’s stated in the error message. Because the current situations is arising when NEV+1 is equal to the global vector length, maybe the check should be updated to accommodate that? Something along the lines of

int needed_vecs = std::min( d_problem->getNEV()+2, MVT::GetGlobalLength(*problem->getInitVec()); TEUCHOS_TEST_FOR_EXCEPTION( needed_vecs > pl.get("Maximum Subspace Dimension”),…);

However you want to modify things to handle this is fine with me.

• Steven

-- Steven Hamilton Radiation Transport Group Reactor and Nuclear Systems Division Oak Ridge National Laboratory Phone: (865) 574-3646_

@aprokop Can you try using Steven's suggested modification to see if it addresses your issues. It should, but I could always be wrong. 😃

[aprokop]

@hkthorn Thank you so much for carefully looking at the issue, and contacting Steven. I'll try to incorporate his suggestion, and will report back.

[mhoemmen]

(listening :) )

[github-actions[bot]]

This issue has had no activity for 365 days and is marked for closure. It will be closed after an additional 30 days of inactivity. If you would like to keep this issue open please add a comment and/or remove the MARKED_FOR_CLOSURE label. If this issue should be kept open even with no activity beyond the time limits you can add the label DO_NOT_AUTOCLOSE. If it is ok for this issue to be closed, feel free to go ahead and close it. Please do not add any comments or change any labels or otherwise touch this issue unless your intention is to reset the inactivity counter for an additional year.

[github-actions[bot]]

This issue was closed due to inactivity for 395 days.