[deal.II] Supported matrix/vector to LinearOperator

[陈敏]

Hi everyone,

When solving a block matrix system, the LinearOperator is a useful tool. My question is whether the LinearOpertor supports Matrix type (the same as the Vector type) like PETSCWrapper::MPI::SparseMatrix or TrilinosWrapper::MPI::SparseMatrix?

Best

Chen

[Bruno Turcksin]

Chen,

Yes, LinearOperator does support Trilinos and PETSc matrices. LinearOperator even supports your own matrix type as long as you define vmult and Tvmult (see here)

Best,

Bruno

[hkch...@gmail.com]

Bruno Turcksin

Thank you! it really helpful!

best

chen

[hkch...@gmail.com]

Hi everyone,

More questions:

1. where can I find the supported linear solver of Trilinos of linear_operator? TrilinosWrappers::SolverDirect seems not to be supported by the linear_operator

2. Does linear_operator can be a preconditioner? for example, SparseMatrix A, B, C, D, E, F; Matrix G=A-B*C^{-1}*D-E^{-1}*F; So Matrix G can only represent by a linear_operator, not its entity. I want to get its inverse by inverse_operator using some linear solvers as: PreconditionILU prec_G(op_G); const auto op_G_inv = inverse_operator(G, Solver, prec_G); But it seems that the precondition of linear_operator is not supported. Is there any solution?

best 

Chen

[Bruno Turcksin]

Chen,

I am not sure if TrilinosWrappers::SolverDirect is supported but for your second point, you can use the linear_operator as a preconditioner. That's exactly what's done here for example. The preconditioner is wrapped in a LinearOperator and this LinearOperator is then used as preconditioner.

Best,

Bruno

[Jean-Paul Pelteret]

Dear Chen,

Bruno pointed you to the exact example that I was going to. It is certainly possible to use linear operators as preconditioners for solvers.

It looks like support for TrilinosWrappers::SolverDirect has not yet been implemented:

https://github.com/dealii/dealii/blob/master/tests/lac/linear_operator_10.cc#L371-L414

I cannot recall what the reason was — maybe the class is simply missing some functions that are defined for the iterative solvers?

As for one of your early questions which Bruno did already answer: I must admit that I do not know if the PETSc linear algebra is compatible with linear operators. We had to write a special TrilinosPayload for the Trilinos linear algebra classes in order for them to work correctly (with the Trilinos solvers, at least), and such a wrapper doesn’t exist for PETSc. But maybe it works just fine if you stick to the deal.II solvers? If you try this out, it would be interesting to know if it works or not.

Best,

J-P

-- 
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 on the web visit https://groups.google.com/d/msgid/dealii/6ee40751-54ae-402d-903e-ae745545387bn%40googlegroups.com.

[hkch...@gmail.com]

Thanks, Bruno and Jean-Paul. I will not use the TrilinosWrapper::SolverDirect in inverse_operator, replaced by an iterative solver. But Bruno, this example initializes a precondition by a matrix 

preconditioner.initialize(system_matrix, data);

While I can't get the entity of the matrix, I want to get the precondition directly from a linear_operator. 

SparseMatrix M; 

auto op_M = linear_operator(M); 

TrilinosWrapper::PreconditionILU prec_M(op_M); 

Or in other hands, Can I get the entity of a matrix from the linear_operator? For example: 

SparseMatrix A, B, C, D, E, F;

linear_operator op_A, op_B, op_C, op_D, op_E, op_F;

inverse_operator  inv_C(op_C, solver, prec_C), inv_E(op_E, solver, prec_E);

linear_operator op_G=op_A-op_B*inv_C*op_D-inv_E*op_F;

Can I get the entity of matrix G by TrilinosWrapper::MPI::SparseMatrix G = op_G; or the precondition of op_G?

auto prec_G = linear_operator(op_G, TrilinosWrapper::PreconditionILU); seems don't wok.

Best,

Chen

[Luca Heltai]

The issue is within Trilinos itself.

The TrilinosPayload (thanks JP!) provides an interface to Epetra_LinearOperator (Trilinos version of our LinearOperator). However, it cannot access matrix elements, as it is not thought for matrix based operations.

You can, however, create a direct solver, wrap it into a LinearOperator, and use the direct solver as a preconditioner, or as one block of a Block preconditioner/solver. Simple workflow:

YourDirectSolverHere Ainv;

Ainv.initialze(A); // A is the actual Trilinos Matrix here, not a linear operator, since direct solvers need access to the matrix entries

auto Aop = linear_operator(A);
auto AinvOp = linear_operator(A, Ainv);

Now you can use Aop and AinvOp, i.e.,

auto S = Bt*AinvOp*B;

would crate a Schur complement.

Take a look at step-70, in the solve() function. There you have an example on how to use them with a block system.

L.

> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/3CCBC467-91F0-47CB-8F49-465FE5DB49E5%40gmail.com.

[Luca Heltai]

> While I can't get the entity of the matrix, I want to get the precondition directly from a linear_operator.


> SparseMatrix M;
> auto op_M = linear_operator(M);
> TrilinosWrapper::PreconditionILU prec_M(op_M);

This you cannot do. It is the other way around. Preconditioners need matrices (linear operator only store information about the *action* of the matrix, or of the various operations. Think of it this way: A*B*C*D_inv is still a linear operator. What matrix should this return to you?)

This is how you should use it, if you want your preconditioner to be a linear operator:

TrilinosWrapper::PreconditionILU prec_M_mat(M);
auto op prec_M = linear_operator(M, prec_M_mat);

> Can I get the entity of matrix G by TrilinosWrapper::MPI::SparseMatrix G = op_G; or the precondition of op_G?

Nope. This, again, does not make sense.

A Matrix is a linear operator (therefore you can create a linear operator from a matrix), but a linear operator is just an action (therefore you cannot create a matrix from a linear operator).

L

[陈敏]

Thanks, Luca.

I learn a lot. Since I use TrilinosWrapper::MPI::Vector, TrilinosWrapper::SolveDirect can't use in linear_operator.

Since I can't get the precondition from linear_operator directly, I will creat a inverse_operator without precondition.

Best,

Chen

Luca Heltai <luca....@gmail.com> 于2022年4月20日周三 15:25写道：

--
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 on the web visit https://groups.google.com/d/msgid/dealii/7126F8FD-80F7-4766-B3C3-985274A281F0%40gmail.com.