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 namespace Thyra {
54 
55  using Teuchos::RCP;
56  using Teuchos::rcp;
57  using Teuchos::ParameterList;
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,...)
94  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
95  typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
96  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
97  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
98  typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
99  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
100  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
101  typedef Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
102  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
103 #ifdef HAVE_MUELU_TPETRA
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();
119  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
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) {
131  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
132  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
133  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
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
148  RCP<XpMat> A00 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
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  }
171  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
172 
173  // Retrieve concrete preconditioner object
174  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
175  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
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
182  RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = Teuchos::null;
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  Teuchos::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) {
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);
204  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
205  }
206  }
207 #endif
208  // build a new MueLu hierarchy
209  H = MueLu::CreateXpetraPreconditioner(A, paramList, coordinates, nullspace);
210 
211  } else {
212  // reuse old MueLu hierarchy stored in MueLu Tpetra/Epetra operator and put in new matrix
213 
214  // get old MueLu hierarchy
215 #if defined(HAVE_MUELU_TPETRA)
216  if (bIsTpetra) {
217 
218  RCP<ThyTpLinOp> tpetr_precOp = rcp_dynamic_cast<ThyTpLinOp>(thyra_precOp);
219  RCP<MueTpOp> muelu_precOp = rcp_dynamic_cast<MueTpOp>(tpetr_precOp->getTpetraOperator(),true);
220 
221  H = muelu_precOp->GetHierarchy();
222  }
223 #endif
224  // TODO add the blocked matrix case here...
225 
226  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
227  "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
228  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
229  "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
230  RCP<MueLu::Level> level0 = H->GetLevel(0);
231  RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
232  RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
233 
234  if (!A0.is_null()) {
235  // If a user provided a "number of equations" argument in a parameter list
236  // during the initial setup, we must honor that settings and reuse it for
237  // all consequent setups.
238  A->SetFixedBlockSize(A0->GetFixedBlockSize());
239  }
240 
241  // set new matrix
242  level0->Set("A", A);
243 
244  H->SetupRe();
245  }
246 
247  // wrap hierarchy H in thyraPrecOp
248  RCP<ThyLinOpBase > thyraPrecOp = Teuchos::null;
249 #if defined(HAVE_MUELU_TPETRA)
250  if (bIsTpetra) {
251  RCP<MueTpOp> muelu_tpetraOp = rcp(new MueTpOp(H));
252  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_tpetraOp));
253  RCP<TpOp> tpOp = Teuchos::rcp_dynamic_cast<TpOp>(muelu_tpetraOp);
254  thyraPrecOp = Thyra::createLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpOp);
255  }
256 #endif
257 
258  if(bIsBlocked) {
259  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::nonnull(thyraPrecOp));
260 
262  //typedef Thyra::XpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyXpLinOp; // unused
263  const RCP<MueXpOp> muelu_xpetraOp = rcp(new MueXpOp(H));
264 
265  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getRangeMap());
266  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getDomainMap());
267 
268  RCP <Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xpOp = Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(muelu_xpetraOp);
269  thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace,xpOp);
270  }
271 
272  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
273 
274  defaultPrec->initializeUnspecified(thyraPrecOp);
275 
276  }
277 
278  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
280  uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
281  TEUCHOS_ASSERT(prec);
282 
283  // Retrieve concrete preconditioner object
284  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
285  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
286 
287  if (fwdOp) {
288  // TODO: Implement properly instead of returning default value
289  *fwdOp = Teuchos::null;
290  }
291 
292  if (supportSolveUse) {
293  // TODO: Implement properly instead of returning default value
294  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
295  }
296 
297  defaultPrec->uninitialize();
298  }
299 
300 
301  // Overridden from ParameterListAcceptor
302  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
304  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
305  paramList_ = paramList;
306  }
307 
308  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
310  return paramList_;
311  }
312 
313  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
315  RCP<ParameterList> savedParamList = paramList_;
316  paramList_ = Teuchos::null;
317  return savedParamList;
318  }
319 
320  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
322  return paramList_;
323  }
324 
325  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
327  static RCP<const ParameterList> validPL;
328 
329  if (Teuchos::is_null(validPL))
330  validPL = rcp(new ParameterList());
331 
332  return validPL;
333  }
334 
335  // Public functions overridden from Teuchos::Describable
336  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
338  return "Thyra::MueLuPreconditionerFactory";
339  }
340 } // namespace Thyra
341 
342 #endif // HAVE_MUELU_STRATIMIKOS
343 
344 #endif // ifdef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
void uninitializePrec(PreconditionerBase< Scalar > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList, Teuchos::RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > coords=Teuchos::null, Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > nullspace=Teuchos::null)
static RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOp, PreconditionerBase< Scalar > *prec, const ESupportSolveUse supportSolveUse) const
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOp) const
Teuchos::RCP< PreconditionerBase< Scalar > > createPrec() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() 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()
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.