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 On Wednesday, August 20, 2025 at 4:11:38 PM UTC-4 [email protected] wrote: > 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. > > > > On 20 Aug 2025, at 19:44, Nathaniel Shaffer <[email protected]> wrote: > > > > 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? 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 > . > > -- 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/3c6af8df-2d99-433a-aa3b-2d8f0bbca452n%40googlegroups.com.
