[Trilinos-Users] [EXTERNAL] Re: Xpetra wrappers around existing Tpetra objects

Sat Jun 4 06:27:28 EDT 2016

Hi Martin,

The "common denominator" (i.e. common superclass) for Tpetra::CrsMatrix 
and the Tpetra-based MueLu operator is:

Tpetra::Operator<>

So you should be able to do sth like:

RCP<Tpetra::Operator<ST,LO,GO,NT> > mueLuPreconditioner = 
MueLu::CreateTpetraPreconditioner(A, paramList);
...
problem->setLeftPrec(mueLuPreconditioner);

Note, a Tpetra::CrsMatrix is a subclass of Tpetra::Operator so this works.

Cheers,
Cihan

On 04/06/16 09:35, Martin Vymazal wrote:
> Hello Andrey and Jonathan,
>
>   the example in the user's guide does not compile. Consider this (I
> only copied the essential pieces from page 22 of
> muelu/doc/UsersGuide/mueluguide.pdf):
>
>   Teuchos::RCP<Tpetra::CrsMatrix<> > A; // A points to type
> Tpetra::CrsMatrix
>
>   Teuchos::RCP<MueLu::TpetraOperator> mueLuPreconditioner =
> MueLu::CreateTpetraPreconditioner(A, paramList); // mueLuPreconditioner
> points to Muelu::TpetraOperator
>
>   // ... create X and B ...
>
>   Belos::LinearProblem<> problem( A, X, B ); // Belos::LinearProblem
> does not see to have default template parameters, so I supplied
>                                              // ScalarType = double, MV
> = Tpetra::MultiVector and OP = Tpetra::CrsMatrix
>                                              // Now Belos::LinearProblem
> thinks that the 'OP' template parameter
>                                              // is Tpetra::CrsMatrix ...
>
>
> problem->setLeftPrec(mueLuPreconditioner);  // ... and expects the same
> type when setting the
>                                              //     preconditioner,
> which is not our case!
>
>
> Not surprisingly, I end up with
>
> error: cannot convert ‘MueLu::TpetraOperator<double, int, int,
> Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial> >*’ to ‘const
> Tpetra::CrsMatrix<double, int, int,
> Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial>, false>*’ in
> initialization
>
>
>   In short, A and mueLuPreconditioner are pointers to different types,
> and Belos::LinearProblem does not like that.
>
>   Can I safely cast the RCP from MueLu::TpetraOperator to
> Tpetra::CrsMatrix? Or should I do something else? Again, at least one
> working example would help a lot ...
>
> Best regards,
>
>   Martin
>
>
> On 2016-06-03 20:11, Andrey Prokopenko wrote:
>> One thing I forgot to mention: you can also provide a tolerance to it
>> like this:
>>
>> double tolerance = 1e-8;
>> H->Iterate(rhs, initial_guess, tolerance);
>>
>> or combine max number of iterations and tolerance like this
>>
>> std::pair<int,double> criteria;
>> criteria.first = max_its;
>> criteria.second = tolerance;
>> H->Iterate(rhs, initial_guess, criteria);
>>
>> -Andrey
>>
>> On 06/03/2016 12:42 PM, Martin Vymazal wrote:
>>> Hi Andrey,
>>>
>>>  yes, that's exactly what I meant. Thank you.
>>>
>>> Best regards,
>>>
>>>  Martin
>>>
>>> On 2016-06-03 18:22, Andrey Prokopenko wrote:
>>>> Hi Martin,
>>>>
>>>> If "MueLu as a solver' you mean a standalone solver, then you just
>>>> need to add the following lines to what Jonathan wrote:
>>>>
>>>> Tpetra::MultiVector<> rhs, initial_guess;
>>>> int max_its = 40;
>>>>
>>>> M->IsPreconditioner(false);
>>>> M->Iterate(rhs, initial_guess, max_its);
>>>>
>>>> -Andrey
>>>>
>>>> On 06/03/2016 10:43 AM, Martin Vymazal wrote:
>>>>> Hello Jonathan,
>>>>>
>>>>>  thank you very much, this helps. I think that's all I need, but I
>>>>> was wondering about one more thing - is there a similarly easy way
>>>>> to use MueLu as a solver?
>>>>>
>>>>> Best regards,
>>>>>
>>>>>  Martin
>>>>>
>>>>>
>>>>> On 2016-06-03 17:28, Jonathan Hu wrote:
>>>>>> Hi Martin,
>>>>>>
>>>>>>     Yes, there are easier ways to use MueLu with an existing
>>>>>> Tpetra::CrsMatrix.  You do not need to convert to Xpetra. Something
>>>>>> like this will work in your case:
>>>>>>
>>>>>> Teuchos::RCP<Tpetra::CrsMatrix<> > A;
>>>>>> ParameterList paramList;
>>>>>> //set the MueLu options
>>>>>> Teuchos::RCP<MueLu::TpetraOperator> M =
>>>>>> MueLu::CreateTpetraPreconditioner(A, paramList);
>>>>>>
>>>>>> M is a Tpetra::Operator and can be used directly with Belos.
>>>>>>
>>>>>> This is described in more detail in the user's guide
>>>>>> (muelu/doc/UsersGuide/mueluguide.pdf) on page 22.  Thank you for the
>>>>>> feedback on the documentation ... I'll see if we can make the
>>>>>> migration path clearer for first-time users.
>>>>>>
>>>>>> Regards,
>>>>>> Jonathan
>>>>>>
>>>>>> Martin Vymazal wrote on 06/03/2016 08:10 AM:
>>>>>>> Hello,
>>>>>>>
>>>>>>>  I have an RCP pointer to Tpetra::CrsMatrix called A and I would
>>>>>>> like to pass this to MueLu Hierarchy object. It seems that the
>>>>>>> only option is to wrap A, the right-hand side vector and the
>>>>>>> vector of unknows in Xpetra RCP pointers. The following compiles
>>>>>>> (I found a similar code in the MueLU tutorial):
>>>>>>>
>>>>>>>   Teuchos::RCP<Xpetra::CrsMatrix<>> exA = Teuchos::rcp(new
>>>>>>> Xpetra::TpetraCrsMatrix<>(A));
>>>>>>>   Teuchos::RCP<Xpetra::CrsMatrixWrap<>> exAWrap =
>>>>>>> Teuchos::rcp(new Xpetra::CrsMatrixWrap<>(exA));
>>>>>>>   RCP<Xpetra::Matrix<>> xpetra_A =
>>>>>>> Teuchos::rcp_dynamic_cast<Xpetra::Matrix<>>(exAWrap);
>>>>>>>
>>>>>>>   MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal> H(xpetra_A);
>>>>>>>   H.SetDefaultVerbLevel(MueLu::Medium);
>>>>>>>
>>>>>>> Could you please confirm that the first line only creates a
>>>>>>> shallow copy of A? Is there any easier/more correct/more elegant
>>>>>>> way of constructing xpetra_A? I'm also stuck at this point,
>>>>>>> because now I need Xpetra-style pointers to the map that was used
>>>>>>> to construct (Tpetra) A in order to construct Xpetra-style
>>>>>>> right-hand side multivector (or pointer to it). How do I get
>>>>>>> Xpetra vectors from Tpetra vectors, given the fact that the map
>>>>>>> of A is a Tpetra map?
>>>>>>>
>>>>>>> Thank you for your help.
>>>>>>>
>>>>>>>  Martin
>>>>>>>
>>>>>>> P. S. MueLu looks interesting, but the documentation (especially
>>>>>>> the Xpetra part) could be a bit more extensive ... It looks like
>>>>>>> there's a lot of attention paid to users who migrate from Epetra
>>>>>>> to Tpetra, but what about those, who already have codebase which
>>>>>>> is exclusively Tpetra-based? Are they supposed to move everything
>>>>>>> to Xpetra in order to be able to use MueLu?
>>>>>>> _______________________________________________
>>>>>>> Trilinos-Users mailing list
>>>>>>> Trilinos-Users at trilinos.org
>>>>>>> https://trilinos.org/mailman/listinfo/trilinos-users
>>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> Trilinos-Users mailing list
>>>>>> Trilinos-Users at trilinos.org
>>>>>> https://trilinos.org/mailman/listinfo/trilinos-users
>>>>> _______________________________________________
>>>>> Trilinos-Users mailing list
>>>>> Trilinos-Users at trilinos.org
>>>>> https://trilinos.org/mailman/listinfo/trilinos-users
>>> _______________________________________________
>>> Trilinos-Users mailing list
>>> Trilinos-Users at trilinos.org
>>> https://trilinos.org/mailman/listinfo/trilinos-users
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at trilinos.org
> https://trilinos.org/mailman/listinfo/trilinos-users