Improve UMFPack vector chooser

UMFPack solver can handle nested matrices due to the use of flat-vector/matrix algorithms, however, the UMFPackVectorChooser does not reflect that (e.g., a range type of a BCRSMatrix<BCRSMatrix<...>> is a BlockVector<BCRSMatrix<...>>). This MR, simple makes the UMFPackVectorChooser more friendly to those nested types.

• Solve nested type BCRSMatrix and FieldVector with recursion of UMFPackVectorChooser

• Also, resolve the UMFPackCreator type with SFINAE on the UMFPackVectorChooser instead of arbitrary rules of the creator (the vector chooser already knows the rules, there is no need to repeat them)

• Allow custom vectors in the UMFPack solver

• Add README

changed the description

marked this merge request as draft

added 1 commit

• e7cf80d1 - Solve nested type BCRSMatrix and FieldVector with recursion

Compare with previous version

added 1 commit

• ad5c398f - Add SFINAE requirements to other specialization & allow custom vectors

Compare with previous version

changed the description

added 1 commit

• a3a16886 - Remove default domain/range types from UMFPack signature

Compare with previous version

changed the description

resolved all threads

changed the description

marked this merge request as ready

added feature label

requested review from @oliver.sander

I don't understand the details of this, but I looked at your changes and didn't see anything obviously wrong with it. Since it doesn't break the existing tests it can be merged IMO. Even better would be if you could extend the test such that your new features get tested as well.

And it needs a changelog entry.

My intention was to use this in arbitrary nested matrices (i.e., BCRSMatrix<BCRSMatrix<...>>), but I just learned that flatMatrixForEach only accepts "well aligned blocks" :( The svg writer had the same issue and the solution was to store the offets of each row/column, but I don't think that's performant here, @patrick.jaap perhaps you have an idea for this?

Hi. flatMatrixForEach should work for BCRSMatrix<BCRSMatrix<...>> as long as there is at least one initialized block entry in any nested BCRSMatrix, and all blocks entries within a nested BCRSMatrix have the same size. Can you be more specific what does not work with UMFPack in your case?

[...] all blocks entries within a nested BCRSMatrix have the same size.

This. It only works in few cases. Imagine each block comes from the discretization of a different sub-domain or equation, then each sub-block may have a completely different pattern.

I think using nested BCRSMatrix types with different sizes can create some strange corner cases. For example, you cannot leave a row/col empty, and computing the flat (outer) dimension of the matrix is complicated. But maybe I am too closed-minded here. If your outer BCRSMatrix is dense, then I guess MultiTypeBlockMatrix would be a better choice in your setting? Otherwise we would need a more general implementation of flatMatrixForEach, which does a more flexible iteration of the flat indices and checks the consistency of the blocking all the time.

I was thinking that if one had the domain and range vectors of the linear map it would be very easy to create the offets from domain/range size.

If your outer BCRSMatrix is dense, then I guess MultiTypeBlockMatrix would be a better choice in your setting?

How would that prevent sub-block having different sizes? In either case, I do want it to be sparse and dynamic because that's the setting of the problem.

Otherwise we would need a more general implementation of flatMatrixForEach, which does a more flexible iteration of the flat indices and checks the consistency of the blocking all the time.

This is not impossible as I mentioned above, it's just kind of expensive because you need two passes through the matrix.

When I implemented the more generic data types in UMFPack I had MultiTypeBlockMatrix with different BCRSMatrix<FieldMatrix<double,n,m>> in mind. Therefore, each BCRSMatrix has constant block sizes, and different block sizes are separated by the MultiTypeBlockMatrix. I still find nesting BCRSMatrix types a bit strange. And as you mentioned, it is also not easy to efficiently put this structure into UMFPack.

I still do not see how MultiTypeBlockMatrix of BCRSMatrix is guaranteed to have same block sizes. But in any case, I think the function name over-promises what it can actually delivers. It's not a critique, I think that is fine to assume regular blocks, I just got confused by the name . I would propose to rename it then.

No, it does not guarantee to have same block sizes. That's why I mentioned BCRSMatrix<FieldMatrix>>.

added 1 commit

• 6236fc58 - Add comment on allocator rebind

Compare with previous version

I have started a system pipeline on this (https://gitlab.dune-project.org/infrastructure/dune-nightly-test/-/pipelines/73185) If it passes, I will merge.

changed milestone to %DUNE 2.10.0

mentioned in commit cac86579

merged

@santiago.ospina I'm not sure this MR is correct. We originally add the specialization of the the creator to detect if a solver would be instantiated over something different from double. Looking at the changes in isValidBlock I don't see how you guarantee this now.

mentioned in merge request !453 (merged)

There are no tests that would trigger your problem.

OK, I now see how it is intended to work... You simply check that the two types are the same and in the case of a non-double field the UMFPackVectorChooser does not specify a type and thus the is_same should return false.

Still a test would be highly appreciated to avoid regressions.

Long story short, I didn't have enough time to solve all the problems and make a proper test.

I wanted to do the test, the problem is that I wanted to use this for nested blocks of BCRS matrices (see description). However, as discussed above, the implementation of flatMatrixForEach is not generic enough to account for different block sizes. Allowing for the type is just half of the story in my opinion... I guess that there are other nested types that would have worked, but that was not my goal for this MR.

In the end, I had to roll out my own implementation¹ as a generic solution for flatMatrixForEach would be quite some work and I needed a quick solution that works for the nested types that I use.

If this MR makes you life too hard for !453 (merged), I would propose to just revert it and rethink this together with a solution for flatMatrixForEach that would properly allow to use this wrapper on nested matrices.

1. Luckily, I didn't look at the actual internal implementation of the umpfpack wrapper before I wrote mine, I just eye witnessed that flatMatrixForEach caused the problem. Therefore a different license was possible.

OK, I now see how it is intended to work... You simply check that the two types are the same and in the case of a non-double field the UMFPackVectorChooser does not specify a type and thus the is_same should return false.

Exactly. Note how the previous version was basically repeating the rules of the vector chooser. This way is easier to add more rules to the chooser if that's needed in the future.

I got it fixed, but I'm not sure that I didn't break your code, as there is not test.

Also there are more solvers using suitesparse. It seems that the implementations slowly start diverging?!

But ... now the code doesn't even allow to instatiate UMFPack<M> if the vector classes are not well defined. At least this makes the creator more cumbersome, see my last changes in !453 (merged) .