MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 namespace Thyra {
54 
55  using Teuchos::RCP;
56  using Teuchos::rcp;
58 
59 
60  // Constructors/initializers/accessors
61 
62  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
64  paramList_(rcp(new ParameterList()))
65  {}
66 
67  // Overridden from PreconditionerFactoryBase
68 
69  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70  bool MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
71  const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
72 
73 #ifdef HAVE_MUELU_TPETRA
74  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp)) return true;
75 #endif
76 
77  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp)) return true;
78 
79  return false;
80  }
81 
82 
83  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
85  return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
86  }
87 
88  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
90  initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
91  using Teuchos::rcp_dynamic_cast;
92 
93  // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
96  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
101  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
102  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
103 #ifdef HAVE_MUELU_TPETRA
105  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpOp;
106  typedef Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyTpLinOp;
107 #endif
108 
109  // Check precondition
110  TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
111  TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
112  TEUCHOS_ASSERT(prec);
113 
114  // Create a copy, as we may remove some things from the list
115  ParameterList paramList = *paramList_;
116 
117  // Retrieve wrapped concrete Xpetra matrix from FwdOp
118  const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
120 
121  // Check whether it is Epetra/Tpetra
122  bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
123  bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
124  bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
125  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
126  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked == false);
127  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked == true);
128 
129  RCP<XpMat> A = Teuchos::null;
130  if(bIsBlocked) {
132  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
134 
135  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
136 
137  Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
138  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
139 
140  RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
141  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
142 
143  // MueLu needs a non-const object as input
144  RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
145  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
146 
147  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
149  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
150 
151  RCP<const XpMap> rowmap00 = A00->getRowMap();
152  RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
153 
154  // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
155  RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
156  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
157 
158  // save blocked matrix
159  A = bMat;
160  } else {
161  RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
162  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
163 
164  // MueLu needs a non-const object as input
165  RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
166  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
167 
168  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
169  A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
170  }
172 
173  // Retrieve concrete preconditioner object
174  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
176 
177  // extract preconditioner operator
178  RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
179  thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
180 
181  // Variable for multigrid hierarchy: either build a new one or reuse the existing hierarchy
183 
184  // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
185  // rebuild preconditioner if startingOver == true
186  // reuse preconditioner if startingOver == false
187  const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");
188 
189  if (startingOver == true) {
190  // extract coordinates from parameter list
191  RCP<XpMultVecDouble> coordinates = Teuchos::null;
193 
194  // TODO check for Xpetra or Thyra vectors?
195  RCP<XpMultVec> nullspace = Teuchos::null;
196 #ifdef HAVE_MUELU_TPETRA
197  if (bIsTpetra) {
198  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
199  RCP<tMV> tpetra_nullspace = Teuchos::null;
200  if (paramList.isType<Teuchos::RCP<tMV> >("Nullspace")) {
201  tpetra_nullspace = paramList.get<RCP<tMV> >("Nullspace");
202  paramList.remove("Nullspace");
203  nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_nullspace);
205  }
206  }
207 #endif
208  // build a new MueLu hierarchy
209  ParameterList& userParamList = paramList.sublist("user data");
210  if(Teuchos::nonnull(coordinates)) {
211  userParamList.set<RCP<XpMultVecDouble> >("Coordinates", coordinates);
212  }
213  if(Teuchos::nonnull(nullspace)) {
214  userParamList.set<RCP<XpMultVec> >("Nullspace", nullspace);
215  }
216  H = MueLu::CreateXpetraPreconditioner(A, paramList);
217 
218  } else {
219  // reuse old MueLu hierarchy stored in MueLu Tpetra/Epetra operator and put in new matrix
220 
221  // get old MueLu hierarchy
222 #if defined(HAVE_MUELU_TPETRA)
223  if (bIsTpetra) {
224 
225  RCP<ThyTpLinOp> tpetr_precOp = rcp_dynamic_cast<ThyTpLinOp>(thyra_precOp);
226  RCP<MueTpOp> muelu_precOp = rcp_dynamic_cast<MueTpOp>(tpetr_precOp->getTpetraOperator(),true);
227 
228  H = muelu_precOp->GetHierarchy();
229  }
230 #endif
231  // TODO add the blocked matrix case here...
232 
234  "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
235  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
236  "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
237  RCP<MueLu::Level> level0 = H->GetLevel(0);
238  RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
239  RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
240 
241  if (!A0.is_null()) {
242  // If a user provided a "number of equations" argument in a parameter list
243  // during the initial setup, we must honor that settings and reuse it for
244  // all consequent setups.
245  A->SetFixedBlockSize(A0->GetFixedBlockSize());
246  }
247 
248  // set new matrix
249  level0->Set("A", A);
250 
251  H->SetupRe();
252  }
253 
254  // wrap hierarchy H in thyraPrecOp
255  RCP<ThyLinOpBase > thyraPrecOp = Teuchos::null;
256 #if defined(HAVE_MUELU_TPETRA)
257  if (bIsTpetra) {
258  RCP<MueTpOp> muelu_tpetraOp = rcp(new MueTpOp(H));
259  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_tpetraOp));
260  RCP<TpOp> tpOp = Teuchos::rcp_dynamic_cast<TpOp>(muelu_tpetraOp);
261  thyraPrecOp = Thyra::createLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpOp);
262  }
263 #endif
264 
265  if(bIsBlocked) {
267 
269  //typedef Thyra::XpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyXpLinOp; // unused
270  const RCP<MueXpOp> muelu_xpetraOp = rcp(new MueXpOp(H));
271 
272  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getRangeMap());
273  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getDomainMap());
274 
276  thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace,xpOp);
277  }
278 
280 
281  defaultPrec->initializeUnspecified(thyraPrecOp);
282 
283  }
284 
285  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
287  uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
288  TEUCHOS_ASSERT(prec);
289 
290  // Retrieve concrete preconditioner object
291  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
293 
294  if (fwdOp) {
295  // TODO: Implement properly instead of returning default value
296  *fwdOp = Teuchos::null;
297  }
298 
299  if (supportSolveUse) {
300  // TODO: Implement properly instead of returning default value
301  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
302  }
303 
304  defaultPrec->uninitialize();
305  }
306 
307 
308  // Overridden from ParameterListAcceptor
309  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
312  paramList_ = paramList;
313  }
314 
315  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
317  return paramList_;
318  }
319 
320  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
322  RCP<ParameterList> savedParamList = paramList_;
323  paramList_ = Teuchos::null;
324  return savedParamList;
325  }
326 
327  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
329  return paramList_;
330  }
331 
332  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
334  static RCP<const ParameterList> validPL;
335 
336  if (Teuchos::is_null(validPL))
337  validPL = rcp(new ParameterList());
338 
339  return validPL;
340  }
341 
342  // Public functions overridden from Teuchos::Describable
343  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
345  return "Thyra::MueLuPreconditionerFactory";
346  }
347 } // namespace Thyra
348 
349 #endif // HAVE_MUELU_STRATIMIKOS
350 
351 #endif // ifdef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
T & get(const std::string &name, T def_value)
static RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool nonnull(const std::shared_ptr< T > &p)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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
bool isParameter(const std::string &name) const
bool remove(std::string const &name, bool throwIfNotExists=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList)
Helper function to create a MueLu preconditioner that can be used by Xpetra.Given an Xpetra::Matrix...
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
Teuchos::RCP< PreconditionerBase< Scalar > > createPrec() const
bool isType(const std::string &name) const
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Exception throws to report errors in the internal logical of the program.
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
bool is_null() const
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.