TridiagonalMatrix: change Assert to AssertThrow

[nfehn]

closes #14609

[peterrum]

Looks good to me. But maybe introducing a new assert type (similar to

dealii/include/deal.II/base/exceptions.h

Line 857 in 6984b3b

DeclException0(ExcNotInitialized);

) would make sense.

My (follow-up) question would be about who is responsible to catch the assert if it thrown during the eigenvalue estimation for PreconditionChebyshev? On the one hand, one could do it in use code (by catching the new assert type in addition to others) or alternatively the assert is catched in SolverCG/PreconditionChebyshev during the eigenvalue estimation and only a single assert is thrown.

/rebuild

[masterleinad]

Looks like some compilers are complaining that info is not initialized. Maybe, just do

types::blas_int info = 0;

[bangerth]

By way of context (just working through this myself): The documentation page at https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/lapack-routines/lapack-least-squares-and-eigenvalue-problem/lapack-least-squares-eigenvalue-problem-driver/symmetric-eigenvalue-problems-lapack-driver/stev.html says that the last argument should actually be called ldz and is an input argument. On the other hand, https://www.smcm.iqfr.csic.es/docs/intel/mkl/mkl_manual/lse/functn_stev.htm provides for the Fortran interface, which I believe we must be using here, where the last argument is an output argument info. It should be zero if operations succeed, and nonzero if they don't -- so an internal error is appropriate.

[bangerth]

I'm ok with this patch if you initialize info as mentioned above.

What isn't quite clear to me is what you want to do with this. My understanding is that if info != 0, then operations actually failed (so ExcInternalError is the right category), but it's not clear to me what this actually implies: Does that mean that the data structures were set up incorrectly (so truly an internal error -- no recovery is possible), or that the matrix has specific properties that one could catch and deal with (in which case one should come up with a different error category)? The documentation states

info INTEGER.

If info = 0, the execution is successful.
If info = -i, the i-th parameter had an illegal value.
If info = i, then the algorithm failed to converge; i elements of e did not converge to zero.

This would suggest that do Assert(..., ExcInternalError()) if info<0 but do AssertThrow(..., ExcSomething()) if info>0 so that one would then try to catch the exception. But I'm unclear what "failed to converge" is supposed to mean.

[bangerth]

Anyone want to take a stab at resolving what we should do here? Do we want this patch, or should we close it? I'm unsure about the semantics of the function called, and don't know how to decide.

[nfehn]

I'm ok with this patch if you initialize info as mentioned above.

I could do this if this is the consensus among the deal.II developers.

What isn't quite clear to me is what you want to do with this.

I gave more motivation (than provided here) in the related issue #14609.

If I remember correctly, I did not want more than catching it. In user code, some high-level parameters like time step size, load factors, boundary conditions, etc. (which might change over time in dynamical problems) might decide about whether a (non)linear system of equations is solvable or not. From a user perspective, I see the motivation to be able to identify "invalid states" and then react to this in the user code. From that perspective, I see an improvement by the suggested change, even though I do not have good answers to all the questions raised here. I think a fundamental question is whether deal.II argues that users have to go into debug mode if the simulation of a dynamic problem fails at some point (while running the same code in all the time steps) or not. @kronbichler You suggested to come up with the present PR (see #14609 (comment)), what do you think?

[kronbichler]

Let me add my take on this: I agree that we should follow the documentation of the respective LAPACK function and with what @bangerth said; the case I think we should address is the one with there is no convergence in the LAPACK algorithm. So I think of something similar to this case we do for the matrix inversion:

dealii/source/lac/lapack_full_matrix.cc

Lines 1751 to 1752 in 81ec864

	Assert(info >= 0, ExcInternalError());
	AssertThrow(info == 0, LACExceptions::ExcSingular());

In fact, I believe that this is a case similar to non-convergence in our iterative solvers: It might be that a user code catches the respective exception and then decides to do something else, e.g. choose another algorithm at run time. This would then support what @bangerth stated, we should introduce or use a specific exception type, to make sure the user code catches this precise exception and not some unrelated exception. Intuitively, for the case info > 0 I would not define a new exception but choose SolverControl::NoConvergence. I admit that we are not using this in a SolverControl context, but the underlying algorithm is really an iterative method of the same kind. But I would also be fine with a new class, that we could insert in the following namespace

dealii/include/deal.II/lac/exceptions.h

Lines 24 to 25 in 81ec864

	namespace LACExceptions
	{

As a follow-up PR, we might want to look over LAPACKFullMatrix and see if there are similar cases that could be covered by this exception type.

[bangerth]

What @kronbichler suggests is pretty much exactly what I had in mind. I'd merge a patch that implements that.