[Core][Future] Add linear operators

[rubenzorrilla]

📝 Description
This PR adds two new linear operator classes. The idea is to enable our linear solvers (and thus strategies) to leverage matrix-free implementations. LinearOperator is the base class defining the API. SparseMatrixLinearOperator is a derived class encapsulating a CSR matrix for the matrix-based implementations (what we had so far).

The idea is that our linear solvers always get a linear operator class instead of a matrix. This statement stems from the idea that any linear operator can be represented as a matrix once a basis is fixed. As a consequence, we will of course require a slight modification in the future linear solvers' API with respect to current one (these changes will be done in an upcoming PR) as we will be always passing a linear operator pointer instead of a reference to the system matrix.

The design is also based on the fact that the linear solver is the one that knows what it needs. Hence, for some of the solvers which accept matrix-free implementations (e.g., CG, GMRES) "any" (by any I mean matrix-free or matrix-based) linear operator can be used. On the contrary, if the linear solver explicitly requires a matrix (e.g., LU decomposition), such matrix is encapsulated and retrieved from a "matrix-based" linear operator, i.e., the SparseMatrixLinearOperator. Note that the query method IsMatrixFree() helps in this regard.

Per regards the type of the matrix stored in the SparseMatrixLinearOperator, we use an std::any to enable storing a reference to the matrix type defined in the TLinearAlgebra template argument (see the linear algebra types #14108). This allows to define the getter in the base class. Furthermore, it also makes the SparseMatrixLinearOperator agnostic to the matrix type or, in other words, any matrix type defining the SpMV and TransposeSpMV methods can work with our SparseMatrixLinearOperator. Note that this does not preclude the possibility to use external matrices (e.g., Epetra) provided that we do a Kratos-like backend to them.

Note that #14108 should be merged first.

[loumalouomega]

No python export?

Okay from my side

[matekelemen]

I don't understand what problem std::any solves here, and my opinion is that using it should have an extremely compelling reason.

Also, why deviate from the linear space model that we have right now? I'd imagine that an implementation of a matrix-free linear space could probably do the same as this proposal does.

[matekelemen]

I don't get this entire design.

First of all, std::any immediately kills the type system. It's unreadable, unsafe, and the compiler is also left completely clueless. This entire system would need painfully detailed documentation just because of this std::any, and I don't think I need to explain that our history with docs is not exactly anything to be proud of.

Secondly, even though GetMatrixImpl returns an std::any, this function is templated on the return type, so the user must somehow know what to expect. Why not just template the whole class on TMatrixType if that's the case? For matrix-free versions, you could just specialize the class on a dummy type (e.g.: MatrixFree).

[rubenzorrilla]

I don't get this entire design.

First of all, std::any immediately kills the type system. It's unreadable, unsafe, and the compiler is also left completely clueless. This entire system would need painfully detailed documentation just because of this std::any, and I don't think I need to explain that our history with docs is not exactly anything to be proud of.

It is something very intrusive that, speaking frankly, an average user shouldn't care about. I note that we have several places in which we do very specific uses of void* in Kratos, and this has never been a problem for the community.

Secondly, even though GetMatrixImpl returns an std::any, this function is templated on the return type, so the user must somehow now what to expect. Why not just template the whole class on TMatrixType if that's the case? For matrix-free versions, you could just specialize the class on a dummy type (e.g.: MatrixFree).

That's exactly what I meant when I said that the linear solver must now how to use the linear operator, and this is what it is (you cannot make a matrix-free LU decomposition for instance). The user must know the type as is the one defining which linear algebra types (or traits) wants to use.

In a preliminary design I did something in the line of your suggestion, having a variadic argument for the TMatrixType. After having discussed it together with @pooyan-dadvand and @RiccardoRossi, we all agreed that this would have been very inconvenient for the users. That's why I decided to go this way.

Having a dummy matrix class that is indeed matrix-free is for me a quick and dirty solution.

[rubenzorrilla]

I don't understand what problem std::any solves here, and my opinion is that using it should have an extremely compelling reason.

Thanks to the std::any any array can be used. In a first design we used a variant, but this would require for anyone trying to work with a new array to add it to the variant in the base class. This can rapidly become messy, on top of the fact that most probably would require includes that you don't really want to include in the core.

Per regards the complexity of the std::any, I presume that very few developers will be modifying the places in which it is used. Indeed, I don't really think we will need to modify it. Note that custom matrix-free linear operators will simply derive from the base class to likely implement the SpMV. In other words, the design is conceived to avoid the user to touch any std::any related stuff. Everything we need is already done in here.

Also, why deviate from the linear space model that we have right now? I'd imagine that an implementation of a matrix-free linear space could probably do the same as this proposal does.

This is much easy to answer. The spaces we have today are an extremely rigid implementation that brought us a lot of nightmares during years. Note that most of the things in there, are now done in the new sparse arrays, so this no longer make sense in the new strategies. Besides, trying to create n-different matrix-free linear operators would lead to n-different space implementations, and consequently to n-different template specializations of everything, something that is OFC no go. This is naively handled by polymorphism with this design.

[rubenzorrilla]

No python export?

Okay from my side

Done. Thanks for pointing it out. To be honest, I completely forgot.

[jcotela]

Looks ok to me overall. I have some minor comments, but I realize that up to a degree it is a matter of taste.

• I would prefer a(n abstract) base class instead of making the matrix version of the linear operator derived from the matrix-free one, since now you are forcing yourself to keep the num rows/cols counters in sync with the matrix size (can the user get the matrix and resize it from the reference? I am not familiar with the upcoming matrix types, but that is something that can be done for ublas matrices for example, and it would desync your data unless the user remembers to call the setters)

• I don't really see an use case for the clear method, as opposed to just deleting the linear operator instance and creating a new one later if needed. For the matrix version I think it is actually a bad idea: the linear operator holds a reference to the matrix, it does not own it, so I don't see why it should get an interface to mess with its memory (you can always get the matrix and clear that if that is what you want to do).

[rubenzorrilla]

Looks ok to me overall. I have some minor comments, but I realize that up to a degree it is a matter of taste.

* I would prefer a(n abstract) base class instead of making the matrix version of the linear operator derived from the matrix-free one, since now you are forcing yourself to keep the num rows/cols counters in sync with the matrix size (can the user get the matrix and resize it from the reference? I am not familiar with the upcoming matrix types, but that is something that can be done for ublas matrices for example, and it would desync your data unless the user remembers to call the setters)

* I don't really see an use case for the clear method, as opposed to just deleting the linear operator instance and creating a new one later if needed. For the matrix version I think it is actually a bad idea: the linear operator holds a reference to the matrix, it does not own it, so I don't see why it should get an interface to mess with its memory (you can always get the matrix and clear that if that is what you want to do).

Thanks for having a look. Indeed, we had a discussion yesterday and some of the comments went in this line. I believe that with the modifications to be added, some of your concerns will disappear. About the Clear, you're probably right, it is a bit misleading.

[rubenzorrilla]

Following our last discussion, current version avoids the std::any and relies on the MatrixType defined in the linear algebra traits. We need the base LinearOperator class to define the complete API to leverage polymorphism so we need to define the GetMatrix in there (which of course throws an error). The return type is the MatrixType in the traits, so the price to pay is that we always need to define the type in the traits, even though you do your own traits for your custom matrix-free linear algebra stuff. This is not that bad (we can eventually create a fake empty matrix for so or just use the one for the standard traits).

The merge-base changed after approval.