Re: [deal.II] Linear operators - TrilinosWrappers

[Juan Felipe Giraldo]

Dear Jean-Paul,

Finally it is working! Thank you for your help and clear explanation!   Now
it is a bit more clear for me how these linear operators work here :).

(Regarding the op_M operator,  I was not using that variable because I
declared M directly as the linear operator of the block(0,0). But thanks to
let me know to avoid problems!)

kind regards,
Felipe


On Wed, Jun 16, 2021 at 7:22 PM Jean-Paul Pelteret 
wrote:

> Dear Felipe,
>
> Changing this
> const auto op_pre = linear_operator(prec_M);
> to this
> const auto op_pre = linear_operator(M, prec_M);
> allowed your program "FEMDG-Dar-par_tofix" to compile for me. If you have
> a look at this documentation for this variant of linear_operator
>
> https://dealii.org/current/doxygen/deal.II/group__LAOperators.html#ga14dbc8c2c27ea3fd45576528a891c6e2
> you'll see that it specifically mentions its use to encapsulate
> preconditioners — this is what I was referring to in my previous message.
> It’s not a particularly easy problem to diagnose, so I’ll give some thought
> on what (if anything) we can do to programatically improve this situation
> on our side. In the mean time, I've added a note to the GitHub issue to
> detail this further.
>
> (BTW. You also probably want to change this
> const auto op_M = linear_operator(M, prec_M);
> to this
> const auto op_M = linear_operator(system_matrix.block(0,
> 0));
> if you are to use it, as op_M then acts like a preconditioner for M, which
> I don’t think is what you want.)
>
> I didn’t try to run the code after checking that it compiled, but I hope
> that this now sorts things out for you. If not, just let us know.
>
> Best,
> Jean-Paul
>
> On 16. Jun 2021, at 12:19, Juan Felipe Giraldo 
> wrote:
>
> Dear Jean-Paul,
>
> Thank you for your reply,
>
> I already set up the Identity preconditioner with a block matrix, but the
> problem remains the same. I am not sure this is the problem because I
> didn't need to set up the identity operator to any matrix in the scenario
> where I could invert. I think my problem is in the Schur preconditioner
> operator before to invert.
> I attach the problematic code and the "working code" I am coding. The only
> difference between them is the line 1251:
>
> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S
> ); //(problematic case)
> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S
> );  //(Case I could invert)
>
> (Please excuse me for the extension of the code, but the solver function
> where I have the issue is too short. In the case you consider it necessary,
> I could code a shorter one)
>
> Thank you so much,
>
> Regards,
> Felipe
>
>
>
> El miércoles, 16 de junio de 2021 a las 15:45:27 UTC+8, Jean-Paul Pelteret
> escribió:
>
>> Dear Felipe,
>>
>> It might be that you need to set up the preconditioner operator with an
>> exemplar matrix (since the identity operator doesn't know the size of the
>> range and domain that it's working on).
>>
>> If that's not the issue then could you please try to reproduce this as a
>> minimal example and share it with us? The program doesn't need to produce
>> any meaningful result, but it would be good if it shows both the working
>> scenario and the problematic one. That way we could investigate further and
>> try to figure out what's going wrong here.
>>
>> Best,
>> Jean-Paul
>>
>> Sent from my mobile device. Please excuse my brevity and any typos.
>>
>> On Wed, 16 Jun 2021, 08:50 Juan Felipe Giraldo, 
>> wrote:
>>
>>> Dear Jean-Paul,
>>>
>>> Thank you for your reply and the complete information you provide me.
>>> Effectively, declaring the preconditioner (out-of-line) was one problem,
>>> but the inverse_operator function is still not matching.
>>>
>>> What I just realized is that when I compute the operator as a
>>> multiplication of a linear operator, an inverse operator and a linear
>>> operator ( for instance, when I compute the Schur complement )
>>>
>>> const auto op_S = op_BT * op_M_inv * op_B;
>>>
>>> this operator becomes to the type
>>> 'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
>>> ,
>>>
>>> but when I compute the operator as a multiplication of other linear
>>> operators directly (for instance, when I compute the Schur preconditioner,
>>> following the procedure in step 20 )
>>>
>>> const auto op_pre = linear_operator(prec_M);
>>> const auto op_aS = op_BT * op_pre * op_B;
>>>
>>> this operator becomes to the type
>>> 'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.
>>>
>>> So, If I use the inverse_operator function with the first operator
>>> (op_S), I can obtain the inverse without any problem,
>>> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S
>>> );
>>>
>>> But if I do it with the second case, it doesn't work because the
>>> functions are not matching.
>>> (const auto op_S_inv = inverse_operator(opa_S, solver_S,
>>> preconditioner_S );
>>>
>>> I am not sure what is happening because all the 

Re: [deal.II] Linear operators - TrilinosWrappers

[Jean-Paul Pelteret]

Dear Felipe,

Changing this
const auto op_pre = linear_operator(prec_M);
to this
const auto op_pre = linear_operator(M, prec_M);
allowed your program "FEMDG-Dar-par_tofix" to compile for me. If you have a 
look at this documentation for this variant of linear_operator
https://dealii.org/current/doxygen/deal.II/group__LAOperators.html#ga14dbc8c2c27ea3fd45576528a891c6e2
 

you'll see that it specifically mentions its use to encapsulate preconditioners 
— this is what I was referring to in my previous message. It’s not a 
particularly easy problem to diagnose, so I’ll give some thought on what (if 
anything) we can do to programatically improve this situation on our side. In 
the mean time, I've added a note to the GitHub issue to detail this further.

(BTW. You also probably want to change this
const auto op_M  = linear_operator(M, prec_M);
to this
const auto op_M = linear_operator(system_matrix.block(0, 0));
if you are to use it, as op_M then acts like a preconditioner for M, which I 
don’t think is what you want.)

I didn’t try to run the code after checking that it compiled, but I hope that 
this now sorts things out for you. If not, just let us know.

Best,
Jean-Paul

> On 16. Jun 2021, at 12:19, Juan Felipe Giraldo  wrote:
> 
> Dear Jean-Paul,
> 
> Thank you for your reply,
> 
> I already set up the Identity preconditioner with a block matrix, but the 
> problem remains the same. I am not sure this is the problem because I didn't 
> need to set up the identity operator to any matrix in the scenario where I 
> could invert. I think my problem is in the Schur preconditioner operator 
> before to invert.
> I attach the problematic code and the "working code" I am coding. The only 
> difference between them is the line 1251:
> 
> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S ); 
> //(problematic case)
> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );  
> //(Case I could invert)
> 
> (Please excuse me for the extension of the code, but the solver function 
> where I have the issue is too short. In the case you consider it necessary, I 
> could code a shorter one)
> 
> Thank you so much,
> 
> Regards,
> Felipe
> 
> 
> 
> El miércoles, 16 de junio de 2021 a las 15:45:27 UTC+8, Jean-Paul Pelteret 
> escribió:
> Dear Felipe,
> 
> It might be that you need to set up the preconditioner operator with an 
> exemplar matrix (since the identity operator doesn't know the size of the 
> range and domain that it's working on).
> 
> If that's not the issue then could you please try to reproduce this as a 
> minimal example and share it with us? The program doesn't need to produce any 
> meaningful result, but it would be good if it shows both the working scenario 
> and the problematic one. That way we could investigate further and try to 
> figure out what's going wrong here. 
> 
> Best,
> Jean-Paul
> 
> Sent from my mobile device. Please excuse my brevity and any typos.
> 
> On Wed, 16 Jun 2021, 08:50 Juan Felipe Giraldo,  > wrote:
> Dear Jean-Paul,
> 
> Thank you for your reply and the complete information you provide me. 
> Effectively, declaring the preconditioner (out-of-line) was one problem, but 
> the inverse_operator function is still not matching. 
> 
> What I just realized is that when I compute the operator as a multiplication 
> of a linear operator, an inverse operator and a linear operator ( for 
> instance, when I compute the Schur complement )
> 
> const auto op_S = op_BT * op_M_inv * op_B;
> 
> this operator becomes to the type 
> 'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
>  ,
> 
> but when I compute the operator as a multiplication of other linear operators 
> directly (for instance, when I compute the Schur preconditioner, following 
> the procedure in step 20 )
> 
> const auto op_pre = linear_operator(prec_M);
> const auto op_aS = op_BT * op_pre * op_B;
> 
> this operator becomes to the type 
> 'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.
> 
> So, If I use the inverse_operator function with the first operator (op_S), I 
> can obtain the inverse without any problem,
> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );
> 
> But if I do it with the second case, it doesn't work because the functions 
> are not matching. 
> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S );
> 
> I am not sure what is happening because all the operators are declared as  
> "LA::MPI::Vector".
> If you have any suggestion, would be greatly appreciated. 
> 
> Thank you so much,
> 
> Regards,
> Felipe Giraldo
> 
> 
> 
> El martes, 15 de junio de 2021 a las 19:16:54 UTC+8, Jean-Paul Pelteret 
> escribió:
> Hi again Feilpe,
> 
> Regarding the lack of documentation, I’ve opened an issue on Github to track 
> this. You can find that here: 

Re: [deal.II] Linear operators - TrilinosWrappers

[Juan Felipe Giraldo]

Dear Jean-Paul,

Thank you for your reply,

I already set up the Identity preconditioner with a block matrix, but the 
problem remains the same. I am not sure this is the problem because I 
didn't need to set up the identity operator to any matrix in the scenario 
where I could invert. I think my problem is in the Schur preconditioner 
operator before to invert.
I attach the problematic code and the "working code" I am coding. The only 
difference between them is the line 1251:

(const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S 
); //(problematic case)
(const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S 
);  //(Case I could invert)

(Please excuse me for the extension of the code, but the solver function 
where I have the issue is too short. In the case you consider it necessary, 
I could code a shorter one)

Thank you so much,

Regards,
Felipe



El miércoles, 16 de junio de 2021 a las 15:45:27 UTC+8, Jean-Paul Pelteret 
escribió:

> Dear Felipe,
>
> It might be that you need to set up the preconditioner operator with an 
> exemplar matrix (since the identity operator doesn't know the size of the 
> range and domain that it's working on).
>
> If that's not the issue then could you please try to reproduce this as a 
> minimal example and share it with us? The program doesn't need to produce 
> any meaningful result, but it would be good if it shows both the working 
> scenario and the problematic one. That way we could investigate further and 
> try to figure out what's going wrong here. 
>
> Best,
> Jean-Paul
>
> Sent from my mobile device. Please excuse my brevity and any typos.
>
> On Wed, 16 Jun 2021, 08:50 Juan Felipe Giraldo,  
> wrote:
>
>> Dear Jean-Paul,
>>
>> Thank you for your reply and the complete information you provide me. 
>> Effectively, declaring the preconditioner (out-of-line) was one problem, 
>> but the inverse_operator function is still not matching. 
>>
>> What I just realized is that when I compute the operator as a 
>> multiplication of a linear operator, an inverse operator and a linear 
>> operator ( for instance, when I compute the Schur complement )
>>
>> const auto op_S = op_BT * op_M_inv * op_B;
>>
>> this operator becomes to the type 
>> 'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
>>  
>> ,
>>
>> but when I compute the operator as a multiplication of other linear 
>> operators directly (for instance, when I compute the Schur preconditioner, 
>> following the procedure in step 20 )
>>
>> const auto op_pre = linear_operator(prec_M);
>> const auto op_aS = op_BT * op_pre * op_B;
>>
>> this operator becomes to the type 
>> 'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.
>>
>> So, If I use the inverse_operator function with the first operator 
>> (op_S), I can obtain the inverse without any problem,
>> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S 
>> );
>>
>> But if I do it with the second case, it doesn't work because the 
>> functions are not matching. 
>> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S 
>> );
>>
>> I am not sure what is happening because all the operators are declared 
>> as  "LA::MPI::Vector".
>> If you have any suggestion, would be greatly appreciated. 
>>
>> Thank you so much,
>>
>> Regards,
>> Felipe Giraldo
>>
>>
>>
>> El martes, 15 de junio de 2021 a las 19:16:54 UTC+8, Jean-Paul Pelteret 
>> escribió:
>>
>>> Hi again Feilpe,
>>>
>>> Regarding the lack of documentation, I’ve opened an issue on Github to 
>>> track this. You can find that here: 
>>> https://github.com/dealii/dealii/issues/12466
>>>
>>> Best,
>>> Jean-Paul
>>>
>>> On 15. Jun 2021, at 12:57, Jean-Paul Pelteret  
>>> wrote:
>>>
>>> Hi Feilpe,
>>>
>>> Firstly, I agree that the documentation is very light on details on how 
>>> to use the linear operators with Trilinos linear algebra types. We can 
>>> definitely improve this, and any contributions to that effect would be 
>>> greatly appreciated!
>>>
>>> Right now, I can direct you to a few of the tests that use Trilinos LA 
>>> in conjunction with the inverse operator, so that you can compare what you 
>>> have and figure out what the problematic differences are. There is this 
>>> one, for instance
>>>
>>> https://github.com/dealii/dealii/blob/master/tests/lac/linear_operator_12a.cc#L323-L341
>>> that looks like a similar setup to yours, and
>>>
>>> https://github.com/dealii/dealii/blob/master/tests/lac/schur_complement_04.cc#L126-L137
>>> that uses a Trilinos::SolverCG (as opposed to deal.II’s solver).
>>>
>>> Comparing to both of these, I think that the important point might be 
>>> that the preconditioner must be declared (out-of-line) before the 
>>> inverse_operation() function is called. The linear operators typically 
>>> expect the lifetime of the LA objects to exceed that of the operators, and 
>>> so you have to first create the matrix or preconditioner and then pass it 
>>> 

Re: [deal.II] Linear operators - TrilinosWrappers

[Jean-Paul Pelteret]

Dear Felipe,

It might be that you need to set up the preconditioner operator with an
exemplar matrix (since the identity operator doesn't know the size of the
range and domain that it's working on).

If that's not the issue then could you please try to reproduce this as a
minimal example and share it with us? The program doesn't need to produce
any meaningful result, but it would be good if it shows both the working
scenario and the problematic one. That way we could investigate further and
try to figure out what's going wrong here.

Best,
Jean-Paul

Sent from my mobile device. Please excuse my brevity and any typos.

On Wed, 16 Jun 2021, 08:50 Juan Felipe Giraldo, 
wrote:

> Dear Jean-Paul,
>
> Thank you for your reply and the complete information you provide me.
> Effectively, declaring the preconditioner (out-of-line) was one problem,
> but the inverse_operator function is still not matching.
>
> What I just realized is that when I compute the operator as a
> multiplication of a linear operator, an inverse operator and a linear
> operator ( for instance, when I compute the Schur complement )
>
> const auto op_S = op_BT * op_M_inv * op_B;
>
> this operator becomes to the type
> 'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
> ,
>
> but when I compute the operator as a multiplication of other linear
> operators directly (for instance, when I compute the Schur preconditioner,
> following the procedure in step 20 )
>
> const auto op_pre = linear_operator(prec_M);
> const auto op_aS = op_BT * op_pre * op_B;
>
> this operator becomes to the type
> 'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.
>
> So, If I use the inverse_operator function with the first operator (op_S),
> I can obtain the inverse without any problem,
> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );
>
> But if I do it with the second case, it doesn't work because the functions
> are not matching.
> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S
> );
>
> I am not sure what is happening because all the operators are declared as
> "LA::MPI::Vector".
> If you have any suggestion, would be greatly appreciated.
>
> Thank you so much,
>
> Regards,
> Felipe Giraldo
>
>
>
> El martes, 15 de junio de 2021 a las 19:16:54 UTC+8, Jean-Paul Pelteret
> escribió:
>
>> Hi again Feilpe,
>>
>> Regarding the lack of documentation, I’ve opened an issue on Github to
>> track this. You can find that here:
>> https://github.com/dealii/dealii/issues/12466
>>
>> Best,
>> Jean-Paul
>>
>> On 15. Jun 2021, at 12:57, Jean-Paul Pelteret  wrote:
>>
>> Hi Feilpe,
>>
>> Firstly, I agree that the documentation is very light on details on how
>> to use the linear operators with Trilinos linear algebra types. We can
>> definitely improve this, and any contributions to that effect would be
>> greatly appreciated!
>>
>> Right now, I can direct you to a few of the tests that use Trilinos LA in
>> conjunction with the inverse operator, so that you can compare what you
>> have and figure out what the problematic differences are. There is this
>> one, for instance
>>
>> https://github.com/dealii/dealii/blob/master/tests/lac/linear_operator_12a.cc#L323-L341
>> that looks like a similar setup to yours, and
>>
>> https://github.com/dealii/dealii/blob/master/tests/lac/schur_complement_04.cc#L126-L137
>> that uses a Trilinos::SolverCG (as opposed to deal.II’s solver).
>>
>> Comparing to both of these, I think that the important point might be
>> that the preconditioner must be declared (out-of-line) before the
>> inverse_operation() function is called. The linear operators typically
>> expect the lifetime of the LA objects to exceed that of the operators, and
>> so you have to first create the matrix or preconditioner and then pass it
>> to the linear operators. This is similar to what you’d done when setting up 
>> op_M,
>> for example. The case of passing in a deal::PreconditionerIdentity() for
>> serial operations is a special case, and I don’t think that we’ve
>> duplicated that for TrilinosWrappers:: PreconditionerIdentity(). Maybe
>> that could be improved too.
>>
>> I hope that with this single change you’d be able to get your
>> program running. If not, then please do let us know so that we can try to
>> help further.
>>
>> Best,
>> Jean-Paul
>>
>>
>> On 14. Jun 2021, at 19:18, Juan Felipe Giraldo 
>> wrote:
>>
>> Hello everyone,
>>
>> I have implemented a residual-minimization framework that somehow is
>> similar to DPG. I want to extend my results by using parallelization using
>> MPI with PETSc or Trilinos.
>> So far, I have solved the saddle point problem using the Schur complement
>> exactly how it is described in step 20. Now, I am trying to replicate
>> exactly the same solver but using the MPI wrappers and the linear operators.
>>
>> The problem is that when I am trying to implement the inverse_operator to
>> compute the Preconditioner of the Schur complement, I 

Re: [deal.II] Linear operators - TrilinosWrappers

[Juan Felipe Giraldo]

Dear Jean-Paul,

Thank you for your reply and the complete information you provide me. 
Effectively, declaring the preconditioner (out-of-line) was one problem, 
but the inverse_operator function is still not matching. 

What I just realized is that when I compute the operator as a 
multiplication of a linear operator, an inverse operator and a linear 
operator ( for instance, when I compute the Schur complement )

const auto op_S = op_BT * op_M_inv * op_B;

this operator becomes to the type 
'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
 
,

but when I compute the operator as a multiplication of other linear 
operators directly (for instance, when I compute the Schur preconditioner, 
following the procedure in step 20 )

const auto op_pre = linear_operator(prec_M);
const auto op_aS = op_BT * op_pre * op_B;

this operator becomes to the type 
'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.

So, If I use the inverse_operator function with the first operator (op_S), 
I can obtain the inverse without any problem,
(const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );

But if I do it with the second case, it doesn't work because the functions 
are not matching. 
(const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S );

I am not sure what is happening because all the operators are declared as  
"LA::MPI::Vector".
If you have any suggestion, would be greatly appreciated. 

Thank you so much,

Regards,
Felipe Giraldo



El martes, 15 de junio de 2021 a las 19:16:54 UTC+8, Jean-Paul Pelteret 
escribió:

> Hi again Feilpe,
>
> Regarding the lack of documentation, I’ve opened an issue on Github to 
> track this. You can find that here: 
> https://github.com/dealii/dealii/issues/12466
>
> Best,
> Jean-Paul
>
> On 15. Jun 2021, at 12:57, Jean-Paul Pelteret  wrote:
>
> Hi Feilpe,
>
> Firstly, I agree that the documentation is very light on details on how to 
> use the linear operators with Trilinos linear algebra types. We can 
> definitely improve this, and any contributions to that effect would be 
> greatly appreciated!
>
> Right now, I can direct you to a few of the tests that use Trilinos LA in 
> conjunction with the inverse operator, so that you can compare what you 
> have and figure out what the problematic differences are. There is this 
> one, for instance
>
> https://github.com/dealii/dealii/blob/master/tests/lac/linear_operator_12a.cc#L323-L341
> that looks like a similar setup to yours, and
>
> https://github.com/dealii/dealii/blob/master/tests/lac/schur_complement_04.cc#L126-L137
> that uses a Trilinos::SolverCG (as opposed to deal.II’s solver).
>
> Comparing to both of these, I think that the important point might be that 
> the preconditioner must be declared (out-of-line) before the 
> inverse_operation() function is called. The linear operators typically 
> expect the lifetime of the LA objects to exceed that of the operators, and 
> so you have to first create the matrix or preconditioner and then pass it 
> to the linear operators. This is similar to what you’d done when setting up 
> op_M, 
> for example. The case of passing in a deal::PreconditionerIdentity() for 
> serial operations is a special case, and I don’t think that we’ve 
> duplicated that for TrilinosWrappers:: PreconditionerIdentity(). Maybe 
> that could be improved too.
>
> I hope that with this single change you’d be able to get your 
> program running. If not, then please do let us know so that we can try to 
> help further.
>
> Best,
> Jean-Paul
>
>
> On 14. Jun 2021, at 19:18, Juan Felipe Giraldo  wrote:
>
> Hello everyone,
>
> I have implemented a residual-minimization framework that somehow is 
> similar to DPG. I want to extend my results by using parallelization using 
> MPI with PETSc or Trilinos. 
> So far, I have solved the saddle point problem using the Schur complement 
> exactly how it is described in step 20. Now, I am trying to replicate 
> exactly the same solver but using the MPI wrappers and the linear operators.
>
> The problem is that when I am trying to implement the inverse_operator to 
> compute the Preconditioner of the Schur complement, I get an error saying 
> that the functions are not matching "inverse_operator(op_aS, solver_aS, 
> TrilinosWrappers::PreconditionIdentity())."
>
> There is no much documentation about linear operators in parallel solvers, 
> so if anyone has any suggestion on how to fix this problem, it would be 
> well appreciated. 
>
> I have pasted the complete function in below:
>
>
> template 
> void FEMwDG::
>   solve()
>   {
> TimerOutput::Scope t(computing_timer, "solve");
>
> LA::MPI::PreconditionJacobi prec_M;
>  LA::MPI::PreconditionJacobi::AdditionalData data;
>  prec_M.initialize(system_matrix.block(0, 0), data);
> }
>
> auto  = solution.block(0);
> auto  = solution.block(1);
> const auto  = system_rhs.block(0);
> const auto  M = 
> 

Re: [deal.II] Linear operators - TrilinosWrappers

[Jean-Paul Pelteret]

Hi again Feilpe,

Regarding the lack of documentation, I’ve opened an issue on Github to track 
this. You can find that here: https://github.com/dealii/dealii/issues/12466 


Best,
Jean-Paul

> On 15. Jun 2021, at 12:57, Jean-Paul Pelteret  wrote:
> 
> Hi Feilpe,
> 
> Firstly, I agree that the documentation is very light on details on how to 
> use the linear operators with Trilinos linear algebra types. We can 
> definitely improve this, and any contributions to that effect would be 
> greatly appreciated!
> 
> Right now, I can direct you to a few of the tests that use Trilinos LA in 
> conjunction with the inverse operator, so that you can compare what you have 
> and figure out what the problematic differences are. There is this one, for 
> instance
> https://github.com/dealii/dealii/blob/master/tests/lac/linear_operator_12a.cc#L323-L341
>  
> 
> that looks like a similar setup to yours, and
> https://github.com/dealii/dealii/blob/master/tests/lac/schur_complement_04.cc#L126-L137
>  
> 
> that uses a Trilinos::SolverCG (as opposed to deal.II’s solver).
> 
> Comparing to both of these, I think that the important point might be that 
> the preconditioner must be declared (out-of-line) before the 
> inverse_operation() function is called. The linear operators typically expect 
> the lifetime of the LA objects to exceed that of the operators, and so you 
> have to first create the matrix or preconditioner and then pass it to the 
> linear operators. This is similar to what you’d done when setting up op_M, 
> for example. The case of passing in a deal::PreconditionerIdentity() for 
> serial operations is a special case, and I don’t think that we’ve duplicated 
> that for TrilinosWrappers:: PreconditionerIdentity(). Maybe that could be 
> improved too.
> 
> I hope that with this single change you’d be able to get your program 
> running. If not, then please do let us know so that we can try to help 
> further.
> 
> Best,
> Jean-Paul
> 
> 
>> On 14. Jun 2021, at 19:18, Juan Felipe Giraldo > > wrote:
>> 
>> Hello everyone,
>> 
>> I have implemented a residual-minimization framework that somehow is similar 
>> to DPG. I want to extend my results by using parallelization using MPI with 
>> PETSc or Trilinos. 
>> So far, I have solved the saddle point problem using the Schur complement 
>> exactly how it is described in step 20. Now, I am trying to replicate 
>> exactly the same solver but using the MPI wrappers and the linear operators.
>> 
>> The problem is that when I am trying to implement the inverse_operator to 
>> compute the Preconditioner of the Schur complement, I get an error saying 
>> that the functions are not matching "inverse_operator(op_aS, solver_aS, 
>> TrilinosWrappers::PreconditionIdentity())."
>> 
>> There is no much documentation about linear operators in parallel solvers, 
>> so if anyone has any suggestion on how to fix this problem, it would be well 
>> appreciated. 
>> 
>> I have pasted the complete function in below:
>> 
>> 
>> template 
>> void FEMwDG::
>>   solve()
>>   {
>> TimerOutput::Scope t(computing_timer, "solve");
>> 
>> LA::MPI::PreconditionJacobi prec_M;
>>  LA::MPI::PreconditionJacobi::AdditionalData data;
>>  prec_M.initialize(system_matrix.block(0, 0), data);
>> }
>> 
>> auto  = solution.block(0);
>> auto  = solution.block(1);
>> const auto  = system_rhs.block(0);
>> const auto  M = linear_operator(system_matrix.block(0, 
>> 0));
>> 
>> const auto op_M  = linear_operator(M, prec_M);
>> const auto op_B   = 
>> linear_operator(system_matrix.block(0, 1));
>> const auto op_BT = 
>> linear_operator(system_matrix.block(1, 0));
>> 
>> ReductionControl  inner_solver_control2(5000,
>> 1e-18 * 
>> system_rhs.l2_norm(),
>> 1.e-2);
>> 
>> SolverCG cg(inner_solver_control2);
>> const auto op_M_inv = inverse_operator(M, cg, prec_M);
>> 
>> const auto op_S = op_BT * op_M_inv * op_B;
>> const auto op_aS = op_BT * linear_operator(prec_M) * 
>> op_B;
>> 
>> IterationNumberControl   iteration_number_control_aS(30, 1.e-18);
>> SolverCG solver_aS(iteration_number_control_aS);
>> 
>> const auto preconditioner_S =
>> inverse_operator(op_aS, solver_aS, 
>> TrilinosWrappers::PreconditionIdentity());
>> const auto schur_rhs = op_BT * op_M_inv * L ;
>> 
>> SolverControlsolver_control_S(2000, 1.e-12);
>> SolverCG solver_S(solver_control_S);
>> 
>> const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S 
>> );
>> 
>> U = op_S_inv * schur_rhs;
>> std::cout << solver_control_S.last_step()
>>   << 

Re: [deal.II] Linear operators - TrilinosWrappers

[Jean-Paul Pelteret]

Hi Feilpe,

Firstly, I agree that the documentation is very light on details on how to use 
the linear operators with Trilinos linear algebra types. We can definitely 
improve this, and any contributions to that effect would be greatly appreciated!

Right now, I can direct you to a few of the tests that use Trilinos LA in 
conjunction with the inverse operator, so that you can compare what you have 
and figure out what the problematic differences are. There is this one, for 
instance
https://github.com/dealii/dealii/blob/master/tests/lac/linear_operator_12a.cc#L323-L341
 

that looks like a similar setup to yours, and
https://github.com/dealii/dealii/blob/master/tests/lac/schur_complement_04.cc#L126-L137
 

that uses a Trilinos::SolverCG (as opposed to deal.II’s solver).

Comparing to both of these, I think that the important point might be that the 
preconditioner must be declared (out-of-line) before the inverse_operation() 
function is called. The linear operators typically expect the lifetime of the 
LA objects to exceed that of the operators, and so you have to first create the 
matrix or preconditioner and then pass it to the linear operators. This is 
similar to what you’d done when setting up op_M, for example. The case of 
passing in a deal::PreconditionerIdentity() for serial operations is a special 
case, and I don’t think that we’ve duplicated that for TrilinosWrappers:: 
PreconditionerIdentity(). Maybe that could be improved too.

I hope that with this single change you’d be able to get your program running. 
If not, then please do let us know so that we can try to help further.

Best,
Jean-Paul


> On 14. Jun 2021, at 19:18, Juan Felipe Giraldo  wrote:
> 
> Hello everyone,
> 
> I have implemented a residual-minimization framework that somehow is similar 
> to DPG. I want to extend my results by using parallelization using MPI with 
> PETSc or Trilinos. 
> So far, I have solved the saddle point problem using the Schur complement 
> exactly how it is described in step 20. Now, I am trying to replicate exactly 
> the same solver but using the MPI wrappers and the linear operators.
> 
> The problem is that when I am trying to implement the inverse_operator to 
> compute the Preconditioner of the Schur complement, I get an error saying 
> that the functions are not matching "inverse_operator(op_aS, solver_aS, 
> TrilinosWrappers::PreconditionIdentity())."
> 
> There is no much documentation about linear operators in parallel solvers, so 
> if anyone has any suggestion on how to fix this problem, it would be well 
> appreciated. 
> 
> I have pasted the complete function in below:
> 
> 
> template 
> void FEMwDG::
>   solve()
>   {
> TimerOutput::Scope t(computing_timer, "solve");
> 
> LA::MPI::PreconditionJacobi prec_M;
>  LA::MPI::PreconditionJacobi::AdditionalData data;
>  prec_M.initialize(system_matrix.block(0, 0), data);
> }
> 
> auto  = solution.block(0);
> auto  = solution.block(1);
> const auto  = system_rhs.block(0);
> const auto  M = linear_operator(system_matrix.block(0, 
> 0));
> 
> const auto op_M  = linear_operator(M, prec_M);
> const auto op_B   = 
> linear_operator(system_matrix.block(0, 1));
> const auto op_BT = 
> linear_operator(system_matrix.block(1, 0));
> 
> ReductionControl  inner_solver_control2(5000,
> 1e-18 * 
> system_rhs.l2_norm(),
> 1.e-2);
> 
> SolverCG cg(inner_solver_control2);
> const auto op_M_inv = inverse_operator(M, cg, prec_M);
> 
> const auto op_S = op_BT * op_M_inv * op_B;
> const auto op_aS = op_BT * linear_operator(prec_M) * 
> op_B;
> 
> IterationNumberControl   iteration_number_control_aS(30, 1.e-18);
> SolverCG solver_aS(iteration_number_control_aS);
> 
> const auto preconditioner_S =
> inverse_operator(op_aS, solver_aS, 
> TrilinosWrappers::PreconditionIdentity());
> const auto schur_rhs = op_BT * op_M_inv * L ;
> 
> SolverControlsolver_control_S(2000, 1.e-12);
> SolverCG solver_S(solver_control_S);
> 
> const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );
> 
> U = op_S_inv * schur_rhs;
> std::cout << solver_control_S.last_step()
>   << " CG Schur complement iterations to obtain convergence."
>   << std::endl;
> E = op_M_inv * (L - op_B * U);
>  }
> 
> Thank you so much, 
> 
> Regards,
> Felipe
> 
> 
> 
> 
> 
> -- 
> 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 

[deal.II] Linear operators - TrilinosWrappers

[Juan Felipe Giraldo]

Hello everyone,

I have implemented a residual-minimization framework that somehow is 
similar to DPG. I want to extend my results by using parallelization using 
MPI with PETSc or Trilinos. 
So far, I have solved the saddle point problem using the Schur complement 
exactly how it is described in step 20. Now, I am trying to replicate 
exactly the same solver but using the MPI wrappers and the linear operators.

The problem is that when I am trying to implement the inverse_operator to 
compute the Preconditioner of the Schur complement, I get an error saying 
that the functions are not matching "inverse_operator(op_aS, solver_aS, 
TrilinosWrappers::PreconditionIdentity())."

There is no much documentation about linear operators in parallel solvers, 
so if anyone has any suggestion on how to fix this problem, it would be 
well appreciated. 

I have pasted the complete function in below:


template 
void FEMwDG::
  solve()
  {
TimerOutput::Scope t(computing_timer, "solve");

LA::MPI::PreconditionJacobi prec_M;
 LA::MPI::PreconditionJacobi::AdditionalData data;
 prec_M.initialize(system_matrix.block(0, 0), data);
}

auto  = solution.block(0);
auto  = solution.block(1);
const auto  = system_rhs.block(0);
const auto  M = linear_operator(system_matrix.block(0, 
0));

const auto op_M  = linear_operator(M, prec_M);
const auto op_B   = 
linear_operator(system_matrix.block(0, 1));
const auto op_BT = 
linear_operator(system_matrix.block(1, 0));

ReductionControl  inner_solver_control2(5000,
1e-18 * 
system_rhs.l2_norm(),
1.e-2);

SolverCG cg(inner_solver_control2);
const auto op_M_inv = inverse_operator(M, cg, prec_M);

const auto op_S = op_BT * op_M_inv * op_B;
const auto op_aS = op_BT * linear_operator(prec_M) * 
op_B;

IterationNumberControl   iteration_number_control_aS(30, 1.e-18);
SolverCG solver_aS(iteration_number_control_aS);

const auto preconditioner_S =
inverse_operator(op_aS, solver_aS, 
TrilinosWrappers::PreconditionIdentity());
const auto schur_rhs = op_BT * op_M_inv * L ;

SolverControlsolver_control_S(2000, 1.e-12);
SolverCG solver_S(solver_control_S);

const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S 
);

U = op_S_inv * schur_rhs;
std::cout << solver_control_S.last_step()
  << " CG Schur complement iterations to obtain convergence."
  << std::endl;
E = op_M_inv * (L - op_B * U);
 }

Thank you so much, 

Regards,
Felipe




-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/40b8db07-5243-4603-825b-9e7850580206n%40googlegroups.com.