MueLu  Version of the Day
MueLu_TekoSmoother_decl.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 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 
48 #ifndef MUELU_TEKOSMOOTHER_DECL_HPP_
49 #define MUELU_TEKOSMOOTHER_DECL_HPP_
50 
51 #ifdef HAVE_MUELU_TEKO
52 
53 #include "Teko_Utilities.hpp"
54 
55 #include "Teko_InverseLibrary.hpp"
56 #include "Teko_InverseFactory.hpp"
57 
58 #include "MueLu_ConfigDefs.hpp"
59 
60 #include <Teuchos_ParameterList.hpp>
61 
62 #include <Xpetra_MapExtractor_fwd.hpp>
63 
65 #include "MueLu_SmootherPrototype.hpp"
69 #include "MueLu_Monitor.hpp"
70 
71 namespace MueLu {
72 
82  template <class Scalar = SmootherPrototype<>::scalar_type,
83  class LocalOrdinal = typename SmootherPrototype<Scalar>::local_ordinal_type,
86  class TekoSmoother : public SmootherPrototype<Scalar,LocalOrdinal,GlobalOrdinal,Node>
87  {
88  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> MapExtractorClass;
89 
90 #undef MUELU_TEKOSMOOTHER_SHORT
91 #include "MueLu_UseShortNames.hpp"
92 
93  public:
94 
96 
97 
100  TekoSmoother() : type_("Teko smoother") {
101  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::TekoSmoother: Teko can only be used with SC=double. For more information refer to the doxygen documentation of TekoSmoother.");
102  };
103 
105  virtual ~TekoSmoother() { }
107 
109 
110  RCP<const ParameterList> GetValidParameterList() const {
111  RCP<ParameterList> validParamList = rcp(new ParameterList());
112  return validParamList;
113  }
114 
115  void DeclareInput(Level &currentLevel) const { }
116 
117  void SetTekoParameters(RCP<ParameterList> tekoParams) { };
119 
121 
122 
125  void Setup(Level &currentLevel) { }
126 
133  void Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero = false) const { }
135 
136  RCP<SmootherPrototype> Copy() const { return Teuchos::null; }
137 
139 
140 
142  std::string description() const {
143  std::ostringstream out;
145  out << "{type = " << type_ << "}";
146  return out.str();
147  }
148 
150  //using MueLu::Describable::describe; // overloading, not hiding
151  void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel = Default) const {
153 
154  if (verbLevel & Parameters0)
155  out0 << "Prec. type: " << type_ << std::endl;
156 
157  if (verbLevel & Debug)
158  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
159  }
160 
162 
163  private:
164 
166  std::string type_;
167  }; // class TekoSmoother
168 
169 
185  template <class GlobalOrdinal,
186  class Node>
187  class TekoSmoother<double,int,GlobalOrdinal,Node> : public SmootherPrototype<double,int,GlobalOrdinal,Node>
188  {
189  typedef int LocalOrdinal;
190  typedef double Scalar;
191  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> MapExtractorClass;
192 
193 #undef MUELU_TEKOSMOOTHER_SHORT
194 #include "MueLu_UseShortNames.hpp"
195 
196  public:
197 
199 
200 
203  TekoSmoother() : type_("Teko smoother"), A_(Teuchos::null), bA_(Teuchos::null), bThyOp_(Teuchos::null), tekoParams_(Teuchos::null), inverseOp_(Teuchos::null) { };
204 
206  virtual ~TekoSmoother() { }
208 
210 
211  RCP<const ParameterList> GetValidParameterList() const {
212  RCP<ParameterList> validParamList = rcp(new ParameterList());
213 
214  validParamList->set< RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A");
215  validParamList->set< std::string > ("Inverse Type", "", "Name of parameter list within 'Teko parameters' containing the Teko smoother parameters.");
216 
217  return validParamList;
218  }
219 
220  void DeclareInput(Level &currentLevel) const {
221  this->Input(currentLevel, "A");
222  }
223 
224  void SetTekoParameters(RCP<ParameterList> tekoParams) { tekoParams_ = tekoParams; };
226 
228 
229 
232  void Setup(Level &currentLevel) {
233  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
234 
235  FactoryMonitor m(*this, "Setup TekoSmoother", currentLevel);
236  if (this->IsSetup() == true)
237  this->GetOStream(Warnings0) << "MueLu::TekoSmoother::Setup(): Setup() has already been called";
238 
239  // extract blocked operator A from current level
240  A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A"); // A needed for extracting map extractors
241  bA_ = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
242  TEUCHOS_TEST_FOR_EXCEPTION(bA_.is_null(), Exceptions::BadCast,
243  "MueLu::TekoSmoother::Build: input matrix A is not of type BlockedCrsMatrix.");
244 
245  bThyOp_ = bA_->getThyraOperator();
246  TEUCHOS_TEST_FOR_EXCEPTION(bThyOp_.is_null(), Exceptions::BadCast,
247  "MueLu::TekoSmoother::Build: Could not extract thyra operator from BlockedCrsMatrix.");
248 
249  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > thyOp = Teuchos::rcp_dynamic_cast<const Thyra::LinearOpBase<Scalar> >(bThyOp_);
250  TEUCHOS_TEST_FOR_EXCEPTION(thyOp.is_null(), Exceptions::BadCast,
251  "MueLu::TekoSmoother::Build: Downcast of Thyra::BlockedLinearOpBase to Teko::LinearOp failed.");
252 
253  // parameter list contains TekoSmoother parameters but does not handle the Teko parameters itself!
254  const ParameterList& pL = Factory::GetParameterList();
255  std::string smootherType = pL.get<std::string>("Inverse Type");
256  TEUCHOS_TEST_FOR_EXCEPTION(smootherType.empty(), Exceptions::RuntimeError,
257  "MueLu::TekoSmoother::Build: You must provide a 'Smoother Type' name that is defined in the 'Teko parameters' sublist.");
258  type_ = smootherType;
259 
260  TEUCHOS_TEST_FOR_EXCEPTION(tekoParams_.is_null(), Exceptions::BadCast,
261  "MueLu::TekoSmoother::Build: No Teko parameters have been set.");
262 
263  Teuchos::RCP<Teko::InverseLibrary> invLib = Teko::InverseLibrary::buildFromParameterList(*tekoParams_);
264  Teuchos::RCP<Teko::InverseFactory> inverse = invLib->getInverseFactory(smootherType);
265 
266  inverseOp_ = Teko::buildInverse(*inverse, thyOp);
267  TEUCHOS_TEST_FOR_EXCEPTION(inverseOp_.is_null(), Exceptions::BadCast,
268  "MueLu::TekoSmoother::Build: Failed to build Teko inverse operator. Probably a problem with the Teko parameters.");
269 
270  this->IsSetup(true);
271  }
272 
279  void Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero = false) const {
280  TEUCHOS_TEST_FOR_EXCEPTION(this->IsSetup() == false, Exceptions::RuntimeError,
281  "MueLu::TekoSmoother::Apply(): Setup() has not been called");
282 
283  Teuchos::RCP<const Teuchos::Comm<int> > comm = X.getMap()->getComm();
284 
285  Teuchos::RCP<const MapExtractor> rgMapExtractor = bA_->getRangeMapExtractor();
286  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgMapExtractor));
287 
288  // copy initial solution vector X to Ptr<Thyra::MultiVectorBase> YY
289 
290  // create a Thyra RHS vector
291  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > thyB = Thyra::createMembers(Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar> >(bThyOp_->productRange()),Teuchos::as<int>(B.getNumVectors()));
292  Teuchos::RCP<Thyra::ProductMultiVectorBase<Scalar> > thyProdB =
293  Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar> >(thyB);
294  TEUCHOS_TEST_FOR_EXCEPTION(thyProdB.is_null(), Exceptions::BadCast,
295  "MueLu::TekoSmoother::Apply: Failed to cast range space to product range space.");
296 
297  // copy RHS vector B to Thyra::MultiVectorBase thyProdB
298  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::updateThyra(Teuchos::rcpFromRef(B), rgMapExtractor, thyProdB);
299 
300  // create a Thyra SOL vector
301  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > thyX = Thyra::createMembers(Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar> >(bThyOp_->productDomain()),Teuchos::as<int>(X.getNumVectors()));
302  Teuchos::RCP<Thyra::ProductMultiVectorBase<Scalar> > thyProdX =
303  Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar> >(thyX);
304  TEUCHOS_TEST_FOR_EXCEPTION(thyProdX.is_null(), Exceptions::BadCast,
305  "MueLu::TekoSmoother::Apply: Failed to cast domain space to product domain space.");
306 
307  // copy RHS vector X to Thyra::MultiVectorBase thyProdX
308  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::updateThyra(Teuchos::rcpFromRef(X), rgMapExtractor, thyProdX);
309 
310  inverseOp_->apply(
311  Thyra::NOTRANS,
312  *thyB, //const MultiVectorBase<Scalar> &X,
313  thyX.ptr(), //const Ptr<MultiVectorBase<Scalar> > &Y,
314  1.0,
315  0.0);
316 
317  // copy back content of Ptr<Thyra::MultiVectorBase> thyX into X
318  Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > XX =
319  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyX, comm);
320 
321  X.update(Teuchos::ScalarTraits<Scalar>::one(), *XX, Teuchos::ScalarTraits<Scalar>::zero());
322  }
324 
325  RCP<SmootherPrototype> Copy() const { return rcp (new TekoSmoother (*this)); }
326 
328 
329 
331  std::string description() const {
332  std::ostringstream out;
334  out << "{type = " << type_ << "}";
335  return out.str();
336  }
337 
339  //using MueLu::Describable::describe; // overloading, not hiding
340  void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel = Default) const {
342 
343  if (verbLevel & Parameters0)
344  out0 << "Prec. type: " << type_ << std::endl;
345 
346  if (verbLevel & Debug)
347  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
348  }
349 
351 
352  private:
353 
355  std::string type_;
356 
358  RCP<Matrix> A_; // < ! internal blocked operator "A" generated by AFact_
359  RCP<BlockedCrsMatrix> bA_;
360  RCP<const Thyra::BlockedLinearOpBase<Scalar> > bThyOp_;
361 
363  RCP<ParameterList> tekoParams_; // < ! parameter list containing Teko parameters. These parameters are not administrated by the factory and not validated.
364 
365  Teko::LinearOp inverseOp_; // < ! Teko inverse operator
366  }; // class TekoSmoother (specialization on SC=double)
367 } // namespace MueLu
368 
369 #define MUELU_TEKOSMOOTHER_SHORT
370 
371 #endif // HAVE_MUELU_TEKO
372 
373 #endif /* MUELU_TEKOSMOOTHER_DECL_HPP_ */
Important warning messages (one line)
Exception indicating invalid cast attempted.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
RCP< SmootherPrototype > Copy() const
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > MapExtractorClass
bool IsSetup() const
Get the state of a smoother prototype.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
RCP< const ParameterList > GetValidParameterList() const
Input.
Timer to be used in factories. Similar to Monitor but with additional timers.
virtual ~TekoSmoother()
Destructor.
Base class for smoother prototypes.
Print additional debugging information.
Namespace for MueLu classes and methods.
std::string description() const
Return a simple one-line description of this object.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Teko smoother.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void SetTekoParameters(RCP< ParameterList > tekoParams)
std::string description() const
Return a simple one-line description of this object.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
std::string type_
smoother type
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
void DeclareInput(Level &currentLevel) const
Input.
virtual const Teuchos::ParameterList & GetParameterList() const
Print class parameters.
RCP< const ParameterList > GetValidParameterList() const
Input.
void Setup(Level &currentLevel)
Setup routine.
void Input(Level &level, const std::string &varName) const
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Teko smoother.
RCP< const Thyra::BlockedLinearOpBase< Scalar > > bThyOp_
Exception throws to report errors in the internal logical of the program.
Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > MapExtractorClass
virtual std::string description() const
Return a simple one-line description of this object.
Interface to block smoothers in Teko package.