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
On Wednesday, August 20, 2025 at 1:14:43 AM UTC-4 [email protected]
wrote:
> Hi Nathaniel,
>
> Maybe a look at LinearOperator
> <https://dealii.org/9.6.0/doxygen/deal.II/classLinearOperator.html>? 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
>
> 在2025年8月20日星期三 UTC+8 05:21:25<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
>>
>
--
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 [email protected].
To view this discussion visit
https://groups.google.com/d/msgid/dealii/852bfe49-26f2-4ce7-9ad5-5c087274f935n%40googlegroups.com.
