Composing matrix-free operators

[Nathaniel Shaffer]

Howdy,

I am implementing a matrix-free shift-and-invert eigensolver for Ax = λBx, where A and B are both symmetric positive-definite MatrixFreeOperators. In the inverse iteration, I need to do a linear solve with the matrix M = A - σB. Since I can't just pass A - σB into, say, SolverCG::solve(), I need to make a new operator for M, call it ShiftedOperator. I'm aware of two strategies to do this, assuming A and B already exist.

1. Derive a new class from MatrixFreeOperators::Base and implement apply_add() and calculate_diagonal() methods.

2. Create a new class, not derived from anything, that just stores const references to A and B and implements a vmult() method in order to meet the requirements of SolverCG::solve().

In both cases, the actual work is being forwarded to the A and B operators, and the new ShiftedOperator just combines those results as needed. Are there strong technical reasons to prefer one approach over the other? If I only need M for this linear solve, is there added value in having the full MatrixFreeOperator functionality? Perhaps I've overlooked something about preconditioning?

Cheers,

Nathaniel

[Chengjiang Yin]

Hi Nathaniel,

Maybe a look at  LinearOperator? This class can wrap linear combinations of the operators as a new linear operator with storage of a few std::function members, lightweight and easy to use.

Regards,

Chengjiang Yin

[Nathaniel Shaffer]

Hi Chenjian,

Thanks for the tip about LinearOperator! That would be very nice if it works, but I was having trouble. I tried the following (apologies for the incomplete snippet, I hope it communicates the gist)

using VectorType = LinearAlgebra::distributed::Vector<double

MatrixFreeOperators::LaplaceOperator<dim, fe_degree, fe_degree+1, 1, VectorType, VectorizedArray<double>> A;

MatrixFreeOperators::MassOperator<dim, fe_degree, fe_degree+1, 1, VectorType, VectorizedArray<double>> B;

...

double shift = 1.0;

LinearOperator<VectorType> op_A = linear_operator(A);

LinearOperator<VectorType> op_B = linear_operator(B);

LinearOperator<VectorType> op_M = op_A - shift*op_B

The error happens inside the call to linear_operator, where there occurs 

linear_operator.h:1221:69: error: non-const lvalue reference to type 'dealii::LinearAlgebra::distributed::Vector<double>' cannot bind to a value of unrelated type 'dealii::Vector<double>'

 1221 |                 [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },

      |                                       

My reading of this is that the lambda wrapping A.vmult() for some reason expects to work with an ordinary gets confused when I give it a distributed vector. 

And at the end of the day, the LinearOperator stuff is much more flexible than what I tried with the wrapper class, but it is accomplishing the same thing, just with lambdas instead of stored references to A and B, right?

Cheers,

Nathaniel

[Luca Heltai]

Dear Nathaniel,

the linear operator function is templated on the vector type, therefore you have to provide it when calling the function:

auto op_A = linear_operator<VectorType>(A);

You can take a look at the dealii/tests/lac/linear_operator*.cc files to see how this can be used.

Best,
Luca.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion visit https://groups.google.com/d/msgid/dealii/852bfe49-26f2-4ce7-9ad5-5c087274f935n%40googlegroups.com.

[Nathaniel Shaffer]

Dear Luca,

Excellent! Thank you for that tip and the pointer to the tests. Indeed, the LinearOperator works exactly as you describe. This will make life much nicer for me going forward! 

Cheers,

Nathaniel

[Alexander]

Hi Nathaniel

While you have seemingly found an answer, I wanted to add that an alternative strategy is to go via the PETSc interface (that also has matrix-free operators) and eventually use SLEPc, which offers, among others, the shift-and-invert spectral transformation (STSINVERT). 

Best wishes

Alexander