MueLu  Version of the Day
Thyra_MueLuPreconditionerFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
47 #define THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
48 
50 
51 #ifdef HAVE_MUELU_STRATIMIKOS
52 
53 // Stratimikos needs Thyra, so we don't need special guards for Thyra here
54 #include "Thyra_DefaultPreconditioner.hpp"
55 #include "Thyra_BlockedLinearOpBase.hpp"
56 #include "Thyra_XpetraLinearOp.hpp"
57 #ifdef HAVE_MUELU_TPETRA
58 #include "Thyra_TpetraLinearOp.hpp"
59 #include "Thyra_TpetraThyraWrappers.hpp"
60 #endif
61 #ifdef HAVE_MUELU_TPETRA
62 #include "Thyra_EpetraLinearOp.hpp"
63 #endif
64 
65 #include "Teuchos_Ptr.hpp"
66 #include "Teuchos_TestForException.hpp"
67 #include "Teuchos_Assert.hpp"
68 #include "Teuchos_Time.hpp"
69 
70 #include <Xpetra_CrsMatrixWrap.hpp>
71 #include <Xpetra_CrsMatrix.hpp>
72 #include <Xpetra_Matrix.hpp>
73 #include <Xpetra_ThyraUtils.hpp>
74 
75 #include <MueLu_Hierarchy.hpp>
77 #include <MueLu_HierarchyHelpers.hpp>
78 #include <MueLu_ParameterListInterpreter.hpp>
79 #include <MueLu_MLParameterListInterpreter.hpp>
80 #include <MueLu_MasterList.hpp>
81 #include <MueLu_XpetraOperator_decl.hpp> // todo fix me
82 #ifdef HAVE_MUELU_TPETRA
83 #include <MueLu_TpetraOperator.hpp>
84 #endif
85 #ifdef HAVE_MUELU_EPETRA
86 #include <MueLu_EpetraOperator.hpp>
87 #endif
88 
89 namespace Thyra {
90 
91  using Teuchos::RCP;
92  using Teuchos::rcp;
93  using Teuchos::ParameterList;
94 
95 
96  // Constructors/initializers/accessors
97 
98  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
100  paramList_(rcp(new ParameterList()))
101  {}
102 
103  // Overridden from PreconditionerFactoryBase
104 
105  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
106  bool MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
107  const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
108 
109 #ifdef HAVE_MUELU_TPETRA
110  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp)) return true;
111 #endif
112 
113 #ifdef HAVE_MUELU_EPETRA
114  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp)) return true;
115 #endif
116 
117  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp)) return true;
118 
119  return false;
120  }
121 
122 
123  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
125  return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
126  }
127 
128  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
130  initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
131  using Teuchos::rcp_dynamic_cast;
132 
133  // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
134  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
135  typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
136  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
137  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
138  typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
139  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
140  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
141  typedef Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
142  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
143 #ifdef HAVE_MUELU_TPETRA
145  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpOp;
146  typedef Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyTpLinOp;
147 #endif
148 #if defined(HAVE_MUELU_EPETRA) and defined(HAVE_MUELU_SERIAL)
149  typedef MueLu::EpetraOperator MueEpOp;
150  typedef Thyra::EpetraLinearOp ThyEpLinOp;
151 #endif
152 
153  //std::cout << "-======---------------------------------" << std::endl;
154  //std::cout << *paramList_ << std::endl;
155  //std::cout << "-======---------------------------------" << std::endl;
156 
157  // Check precondition
158  TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
159  TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
160  TEUCHOS_ASSERT(prec);
161 
162  // Create a copy, as we may remove some things from the list
163  ParameterList paramList = *paramList_;
164 
165  // Retrieve wrapped concrete Xpetra matrix from FwdOp
166  const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
167  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
168 
169  // Check whether it is Epetra/Tpetra
170  bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
171  bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
172  bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
173  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
174  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked == false);
175  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked == true);
176 
177  RCP<XpMat> A = Teuchos::null;
178  if(bIsBlocked) {
179  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
180  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
181  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
182 
183  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
184 
185  Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
186  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
187 
188  RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
189  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
190 
191  // MueLu needs a non-const object as input
192  RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
193  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
194 
195  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
196  RCP<XpMat> A00 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
197  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
198 
199  RCP<const XpMap> rowmap00 = A00->getRowMap();
200  RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
201 
202  // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
203  RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
204  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
205 
206  // save blocked matrix
207  A = bMat;
208  } else {
209  RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
210  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
211 
212  // MueLu needs a non-const object as input
213  RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
214  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
215 
216  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
217  A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
218  }
219  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
220 
221  // Retrieve concrete preconditioner object
222  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
223  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
224 
225  // extract preconditioner operator
226  RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
227  thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
228 
229  // Variable for multigrid hierarchy: either build a new one or reuse the existing hierarchy
230  RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = Teuchos::null;
231 
232  // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
233  // rebuild preconditioner if startingOver == true
234  // reuse preconditioner if startingOver == false
235  const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");
236 
237  if (startingOver == true) {
238  // extract coordinates from parameter list
239  Teuchos::RCP<XpMultVecDouble> coordinates = Teuchos::null;
240 #ifdef HAVE_MUELU_TPETRA
241  if (bIsTpetra) {
242  // Tpetra does not instantiate on Scalar=float by default, so we must check for this
243  // FIXME This will still break if LO != int or GO != int
244  # if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_MUELU_INST_FLOAT_INT_INT)
245  typedef Tpetra::MultiVector<float, LocalOrdinal, GlobalOrdinal, Node> tfMV;
246  RCP<tfMV> floatCoords = Teuchos::null;
247  # endif
248  typedef Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node> tdMV;
249  RCP<tdMV> doubleCoords = Teuchos::null;
250  if (paramList.isType<RCP<tdMV> >("Coordinates")) {
251  doubleCoords = paramList.get<RCP<tdMV> >("Coordinates");
252  paramList.remove("Coordinates");
253  }
254  # if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_MUELU_INST_FLOAT_INT_INT)
255  else if (paramList.isType<RCP<tfMV> >("Coordinates")) {
256  floatCoords = paramList.get<RCP<tfMV> >("Coordinates");
257  paramList.remove("Coordinates");
258  doubleCoords = rcp(new tdMV(floatCoords->getMap(), floatCoords->getNumVectors()));
259  deep_copy(*doubleCoords, *floatCoords);
260  }
261  # endif
262  if(doubleCoords != Teuchos::null) {
263  coordinates = MueLu::TpetraMultiVector_To_XpetraMultiVector<double,LocalOrdinal,GlobalOrdinal,Node>(doubleCoords);
264  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(coordinates));
265  TEUCHOS_TEST_FOR_EXCEPT(doubleCoords->getNumVectors() != coordinates->getNumVectors());
266  }
267  }
268 #endif // HAVE_MUELU_TPETRA
269 
270 #ifdef HAVE_MUELU_EPETRA
271  if (bIsEpetra) {
272  RCP<Epetra_MultiVector> doubleCoords;
273  if (paramList.isType<RCP<Epetra_MultiVector> >("Coordinates")) {
274  doubleCoords = paramList.get<RCP<Epetra_MultiVector> >("Coordinates");
275  paramList.remove("Coordinates");
276  RCP<Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> > epCoordinates = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node>(doubleCoords));
277  RCP<Xpetra::MultiVector<double,int,int,Node> > epCoordinatesMult = rcp_dynamic_cast<Xpetra::MultiVector<double,int,int,Node> >(epCoordinates);
278  coordinates = rcp_dynamic_cast<Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node> >(epCoordinatesMult);
279  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(coordinates));
280  TEUCHOS_TEST_FOR_EXCEPT(doubleCoords->NumVectors() != Teuchos::as<int>(coordinates->getNumVectors()));
281  }
282  }
283 #endif // HAVE_MUELU_EPETRA
284 
285  // TODO check for Xpetra or Thyra vectors?
286 
287  RCP<XpMultVec> nullspace = Teuchos::null;
288 #ifdef HAVE_MUELU_TPETRA
289  if (bIsTpetra) {
290  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
291  RCP<tMV> tpetra_nullspace = Teuchos::null;
292  if (paramList.isType<Teuchos::RCP<tMV> >("Nullspace")) {
293  tpetra_nullspace = paramList.get<RCP<tMV> >("Nullspace");
294  paramList.remove("Nullspace");
295  nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_nullspace);
296  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
297  }
298  }
299 #endif
300 #ifdef HAVE_MUELU_EPETRA
301  if (bIsEpetra) {
302  RCP<Epetra_MultiVector> epetra_nullspace = Teuchos::null;
303  if (paramList.isType<RCP<Epetra_MultiVector> >("Nullspace")) {
304  epetra_nullspace = paramList.get<RCP<Epetra_MultiVector> >("Nullspace");
305  paramList.remove("Nullspace");
306  RCP<Xpetra::EpetraMultiVectorT<int,Node> > xpEpNullspace = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<int,Node>(epetra_nullspace));
307  RCP<Xpetra::MultiVector<double,int,int,Node> > xpEpNullspaceMult = rcp_dynamic_cast<Xpetra::MultiVector<double,int,int,Node> >(xpEpNullspace);
308  nullspace = rcp_dynamic_cast<XpMultVec>(xpEpNullspaceMult);
309  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
310  }
311  }
312 #endif
313 
314  // build a new MueLu hierarchy
315  H = CreateXpetraPreconditioner(A, paramList, coordinates, nullspace);
316 
317  } else {
318  // reuse old MueLu hierarchy stored in MueLu Tpetra/Epetra operator and put in new matrix
319 
320  // get old MueLu hierarchy
321 #ifdef HAVE_MUELU_TPETRA
322  if (bIsTpetra) {
323 
324  RCP<ThyTpLinOp> tpetr_precOp = rcp_dynamic_cast<ThyTpLinOp>(thyra_precOp);
325  RCP<MueTpOp> muelu_precOp = rcp_dynamic_cast<MueTpOp>(tpetr_precOp->getTpetraOperator(),true);
326 
327  H = muelu_precOp->GetHierarchy();
328  }
329 #endif
330 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_SERIAL)
331  if (bIsEpetra) {
332  RCP<ThyEpLinOp> epetr_precOp = rcp_dynamic_cast<ThyEpLinOp>(thyra_precOp);
333  RCP<MueEpOp> muelu_precOp = rcp_dynamic_cast<MueEpOp>(epetr_precOp->epetra_op(),true);
334 
335  H = rcp_dynamic_cast<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(muelu_precOp->GetHierarchy());
336  }
337 #endif
338  // TODO add the blocked matrix case here...
339 
340  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
341  "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
342  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
343  "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
344  RCP<MueLu::Level> level0 = H->GetLevel(0);
345  RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
346  RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
347 
348  if (!A0.is_null()) {
349  // If a user provided a "number of equations" argument in a parameter list
350  // during the initial setup, we must honor that settings and reuse it for
351  // all consequent setups.
352  A->SetFixedBlockSize(A0->GetFixedBlockSize());
353  }
354 
355  // set new matrix
356  level0->Set("A", A);
357 
358  H->SetupRe();
359  }
360 
361  // wrap hierarchy H in thyraPrecOp
362  RCP<ThyLinOpBase > thyraPrecOp = Teuchos::null;
363 #ifdef HAVE_MUELU_TPETRA
364  if (bIsTpetra) {
365  RCP<MueTpOp> muelu_tpetraOp = rcp(new MueTpOp(H));
366  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_tpetraOp));
367  thyraPrecOp = Thyra::createLinearOp(RCP<TpOp>(muelu_tpetraOp));
368  }
369 #endif
370 
371 #if defined(HAVE_MUELU_EPETRA)
372  if (bIsEpetra) {
373  RCP<MueLu::Hierarchy<double,int,int,Xpetra::EpetraNode> > epetraH =
375  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(epetraH), MueLu::Exceptions::RuntimeError,
376  "Thyra::MueLuPreconditionerFactory: Failed to cast Hierarchy to Hierarchy<double,int,int,Xpetra::EpetraNode>. Epetra runs only on the Serial node.");
377  RCP<MueEpOp> muelu_epetraOp = rcp(new MueEpOp(epetraH));
378  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_epetraOp));
379  // attach fwdOp to muelu_epetraOp to guarantee that it will not go away
380  set_extra_data(fwdOp,"IFPF::fwdOp", Teuchos::inOutArg(muelu_epetraOp), Teuchos::POST_DESTROY,false);
381  RCP<ThyEpLinOp> thyra_epetraOp = Thyra::nonconstEpetraLinearOp(muelu_epetraOp, NOTRANS, EPETRA_OP_APPLY_APPLY_INVERSE, EPETRA_OP_ADJOINT_UNSUPPORTED);
382  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyra_epetraOp));
383  thyraPrecOp = rcp_dynamic_cast<ThyLinOpBase>(thyra_epetraOp);
384  }
385 #endif
386 
387  if(bIsBlocked) {
388  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::nonnull(thyraPrecOp));
389 
391  //typedef Thyra::XpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyXpLinOp; // unused
392  const RCP<MueXpOp> muelu_xpetraOp = rcp(new MueXpOp(H));
393 
394  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getRangeMap());
395  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getDomainMap());
396 
397  RCP <Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xpOp = Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(muelu_xpetraOp);
398  thyraPrecOp = Thyra::xpetraLinearOp(thyraRangeSpace, thyraDomainSpace,xpOp);
399  }
400 
401  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
402 
403  defaultPrec->initializeUnspecified(thyraPrecOp);
404 
405  }
406 
407  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
408  Teuchos::RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
409  CreateXpetraPreconditioner(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > op, const Teuchos::ParameterList& inParamList, Teuchos::RCP<Xpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node> > coords, Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > nullspace) const {
410 
413 
414  Teuchos::ParameterList paramList = inParamList;
415 
416  bool hasParamList = paramList.numParams();
417 
418  RCP<HierarchyManager> mueLuFactory;
419 
420  std::string syntaxStr = "parameterlist: syntax";
421  if (hasParamList && paramList.isParameter(syntaxStr) && paramList.get<std::string>(syntaxStr) == "ml") {
422  paramList.remove(syntaxStr);
424  } else {
425  mueLuFactory = rcp(new MueLu::ParameterListInterpreter <Scalar,LocalOrdinal,GlobalOrdinal,Node>(paramList,op->getDomainMap()->getComm()));
426  }
427 
428  RCP<Hierarchy> H = mueLuFactory->CreateHierarchy();
429  H->setlib(op->getDomainMap()->lib());
430 
431 
432  // Set fine level operator
433  H->GetLevel(0)->Set("A", op);
434 
435  // Set coordinates if available
436  if (coords != Teuchos::null) {
437  H->GetLevel(0)->Set("Coordinates", coords);
438  }
439 
440  // Wrap nullspace if available, otherwise use constants
441  if (nullspace == Teuchos::null) {
442  int nPDE = MueLu::MasterList::getDefault<int>("number of equations");
443  if (paramList.isSublist("Matrix")) {
444  // Factory style parameter list
445  const Teuchos::ParameterList& operatorList = paramList.sublist("Matrix");
446  if (operatorList.isParameter("PDE equations"))
447  nPDE = operatorList.get<int>("PDE equations");
448 
449  } else if (paramList.isParameter("number of equations")) {
450  // Easy style parameter list
451  nPDE = paramList.get<int>("number of equations");
452  }
453 
454  nullspace = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(op->getDomainMap(), nPDE);
455  if (nPDE == 1) {
456  nullspace->putScalar(Teuchos::ScalarTraits<Scalar>::one());
457 
458  } else {
459  for (int i = 0; i < nPDE; i++) {
460  Teuchos::ArrayRCP<Scalar> nsData = nullspace->getDataNonConst(i);
461  for (int j = 0; j < nsData.size(); j++) {
462  GlobalOrdinal GID = op->getDomainMap()->getGlobalElement(j) - op->getDomainMap()->getIndexBase();
463 
464  if ((GID-i) % nPDE == 0)
465  nsData[j] = Teuchos::ScalarTraits<Scalar>::one();
466  }
467  }
468  }
469  }
470  H->GetLevel(0)->Set("Nullspace", nullspace);
471 
472 
473  Teuchos::ParameterList nonSerialList,dummyList;
474  MueLu::ExtractNonSerializableData(paramList, dummyList, nonSerialList);
476 
477  mueLuFactory->SetupHierarchy(*H);
478 
479  return H;
480  }
481 
482  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
484  uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
485  TEUCHOS_ASSERT(prec);
486 
487  // Retrieve concrete preconditioner object
488  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
489  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
490 
491  if (fwdOp) {
492  // TODO: Implement properly instead of returning default value
493  *fwdOp = Teuchos::null;
494  }
495 
496  if (supportSolveUse) {
497  // TODO: Implement properly instead of returning default value
498  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
499  }
500 
501  defaultPrec->uninitialize();
502  }
503 
504 
505  // Overridden from ParameterListAcceptor
506  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
508  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
509  paramList_ = paramList;
510  }
511 
512  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
514  return paramList_;
515  }
516 
517  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
519  RCP<ParameterList> savedParamList = paramList_;
520  paramList_ = Teuchos::null;
521  return savedParamList;
522  }
523 
524  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
526  return paramList_;
527  }
528 
529  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
531  static RCP<const ParameterList> validPL;
532 
533  if (Teuchos::is_null(validPL))
534  validPL = rcp(new ParameterList());
535 
536  return validPL;
537  }
538 
539 
540  // Public functions overridden from Teuchos::Describable
541  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
543  return "Thyra::MueLuPreconditionerFactory";
544  }
545 } // namespace Thyra
546 
547 #endif // HAVE_MUELU_STRATIMIKOS
548 
549 #endif // ifdef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
static void AddNonSerializableDataToHierarchy(HierarchyManager &HM, Hierarchy &H, const ParameterList &nonSerialList)
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOp, PreconditionerBase< Scalar > *prec, const ESupportSolveUse supportSolveUse) const
void uninitializePrec(PreconditionerBase< Scalar > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< Teuchos::ParameterList > paramList_
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOp) const
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Class that accepts ML-style parameters and builds a MueLu preconditioner. This interpreter uses the s...
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &paramList, Teuchos::RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > coords, Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > nullspace) const
Teuchos::RCP< PreconditionerBase< Scalar > > createPrec() const
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
Exception throws to report errors in the internal logical of the program.
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
long ExtractNonSerializableData(const Teuchos::ParameterList &inList, Teuchos::ParameterList &serialList, Teuchos::ParameterList &nonSerialList)
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.