Dealii: Incompatible preconditioners for GMG coarse solvers

[appiazzolla]

Hi All,

The constructor of class MGCoarseGridLACIteration<SolverType, VectorType> instantiates an object of class PointerMatrix<PreconditionerType, VectorType> by providing a reference to a PreconditionerType to the constructor.

Following, members void PointerMatrix<MatrixType, VectorType>::vmult_add(VectorType &dst,const VectorType &src) const and void PointerMatrix<MatrixType, VectorType>::Tvmult_add(VectorType &dst,const VectorType &src) const instantiate void PreconditionerType::vmult_add (VectorType &, const VectorType &) const and void PreconditionerType::Tvmult_add (VectorType &, const VectorType &) const respectively.

Several preconditioners, however, have no member void vmult_add (VectorType &, const VectorType &) const and void Tvmult_add (VectorType &, const VectorType &) const (e.g. PreconditionJacobi< MatrixType >, TrilinosWrappers::PreconditionJacobi, PETScWrappers::PreconditionJacobi).

The compiler complains about the lack of these members. My excuses if this is intended.

[tjhei]

The compiler complains about the lack of these members. My excuses if this is intended.

Yes, I have found the same issue before and no, this is not intended. We should investigate if these function calls are really needed and if yes, if we can add them to the TrilinosWrappers.

[bangerth]

Correct. These functions aren't there because nobody found a need so far. But we'd be quite happy to take a patch! What classes concretely do you have in mind?

[appiazzolla]

Well, as I see it the main problem is users not being able to compile preconditioned GMG iterative coarse solvers today.
Just emptying the body of void PointerMatrix<MatrixType, VectorType>::vmult_add(VectorType &dst,const VectorType &src) const and void PreconditionerType::vmult_add (VectorType &, const VectorType &) const and deleting the names of the arguments works for me...

[appiazzolla]

On the long term I agree, either a class PointerPreconditioner<MatrixType, VectorType> without members vmult_add and Tvmult_add or an implementation of vmult_add and Tvmult_add for _all_ preconditioners could be applicable.

[bangerth]

Well, we've always assumed that all linear operators have a vmult_add function. If there are some that don't, then the right thing to do is to add this function, not change other places. Your solution really just papers over the problem because every caller of PointerMatrix::vmult_add will now get a wrong result.

[appiazzolla]

Well, we've always assumed that all linear operators have a vmult_add function. If there are some that don't, then the right thing to do is to add this function, not change other places.

Agreed.

Your solution really just papers over the problem

Correct.

every caller of PointerMatrix::vmult_add will now get a wrong result.

Unless using ExcNotImplemented()

[bangerth]

On 09/09/2016 04:35 PM, appiazzolla wrote:

every caller of PointerMatrix::vmult_add will now get a wrong result.

Unless using |ExcNotImplemented()|

That's a wrong result as well if you really wanted to have the action computed :-)

[appiazzolla]

Maybe MGCoarseGridLACIteration could use a PointerMatrixAux instead of PointerMatrix.

[tjhei]

Yes, that will work. I think we can make MGCoarseGridLACIteration more flexible. I will prepare a pull request...

[tjhei]

can you take a look at #3109 and tell me what you think? (see the test for an example). I didn't want to always use PointerMatrixAux because it is less efficient if it is not needed. I guess the interface is not ideal, but it works at least!

[appiazzolla]

Works like a charm for me!

[appiazzolla]

Maybe it would be a good idea to add in the documentation that a user trying to use a PreconditionerType which does not implement vmult_add and Tvmult_add should "encapsulate" his preconditioner with a PointerMatrixAux before providing it to MGCoarseGridLACIteration.

[appiazzolla]

Couldn't we just get rid of PointerMatrix and all the PointerSomething family?

Something like this:
mg_coarse.txt

In my implementation template argument deduction is not working but I guess it is possible to do so and even if it is not, I would say this is more generic.

Why does the PointerSomething family exist at all?

[tjhei]

Why does the PointerSomething family exist at all?

I think the reason for it was to hide the concrete types of your matrix and preconditioner. With your change, they become template arguments of the class. I don't see this as a big issue though, so maybe we should do what you suggest. Give me some time to think about it.

[appiazzolla]

I may be horribly wrong but all these MGCoarseGridBase< VectorType > classes look to me just like functors, maybe Multigrid< VectorType > should just expect an std::function<void(const unsigned int,VectorType &,const VectorType &)> as coarse argument, then I could just provide a lambda.

[masterleinad]

@appiazzolla Do you have something in mind you want to do what is not already there?

[kronbichler]

For re-design of what we expect from multigrid algorithms, see also #2534.

[masterleinad]

@kronbichler I guess the interface of MGCoarseGridBase is really minimal and we don't want to change it in

2534, isn't it?

In any case, I am much in favor to not force the user to use some intermediate object if the Preconditioner for the coarse grid does not have a vmult_add.

[appiazzolla]

@masterleinad I was thinking about #2534, multigrid classes have several functor implementations which are not compatible between each other, I see lambdas everywhere.
This issue actually solves a problem @davydden comments about in #2534.

[kronbichler]

@masterleinad I didn't think of changing that, I agree. The point is just that a somewhat unified interface to the multigrid algorithms helps. If we go for templates in one place, we might just as well go with templates for both the smoother, the matrix, and the coarse grid solver. (But I would definitely keep base classes with virtual inheritance at least on the smoother side because they allow to use different smoothers on different levels.) So my post was merely to remind ourselves of the big picture. (I know I'm slow on getting my work for #2534 in... :-)).

[davydden]

@masterleinad :

I guess the interface of MGCoarseGridBase is really minimal and we don't want to change it in

2534, isn't it?

i was suggesting that it's enough to have something with vmult https://github.com/dealii/dealii/issues/2534#issuecomment-216890476 without MGCoarseGridBase.

[appiazzolla]

@masterleinad @davydden I would go further to say you just need an std::function<void(const unsigned int,VectorType &,const VectorType &)>

[tjhei]

In any case, I am much in favor to not force the user to use some intermediate object if the Preconditioner for the coarse grid does not have a vmult_add.

Agreed, we can do this easily without changing the base class because we don't need this PointerMatrix stuff.

The point is just that a somewhat unified interface to the multigrid algorithms helps. If we go for templates in one place, we might just as well go with templates for both the smoother, the matrix, and the coarse grid solver. (But I would definitely keep base classes with virtual inheritance at least on the smoother side because they allow to use different smoothers on different levels.)

So the choice is between a) std::function b) base class (like it is now) and c) template argument. I think you are saying that smoothers should do b), so for consistency I think we should stick with b) for all cases then. I kind of like b) best anyways as long as the interface is minimal: this makes MultiGrid have fewer template parameters and has the interface explicitly specified in code (fewer error messages when you get it wrong). I don't consider a) a good solution, but you can always implement a wrapper class for LinearOperator/std::function on top of b, right?

[bangerth]

Just for the record, a virtual function call will cost a tiny fraction of what it costs to do one matrix-vector product. Efficiency should not be a reason to go with templates over base classes, though there may of course be other considerations.

[appiazzolla]

IMHO keeping a base class and have a fully templated MGCoarseGridLACIteration (maybe with the templates arranged such that argument deduction works) solves this issue and goes in line with the present Multigrid<VectorType> implementation.

In the long term the right answer maybe different but that's what #2534 is about, isn't it?

[kronbichler]

I completely agree with @tjhei regarding the design with base classes. Then, we could, as @appiazzolla suggests, build a fully templated class on top of that. (I do not see how to use template deduction, though, but that is not really a problem IMHO.) I'm also fine with the solution sketched in #3109, so I let @tjhei and @appiazzolla decide what fits better.

[appiazzolla]

A final note on template argument deduction: it seems template argument deduction for constructors will be available in C++17.

[bangerth]

I think you forgot to provide the URL for your link. Do you still have it?

[kronbichler]

I guess it is this one:
http://en.cppreference.com/w/cpp/language/class_template_deduction

[appiazzolla]

@bangerth Sorry! I edited the comment.
http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2016/p0091r2.html

[masterleinad]

Fixed by #3109