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

[Martin Vymazal]

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