MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
61 
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,
84  class GlobalOrdinal = typename SmootherPrototype<Scalar, LocalOrdinal>::global_ordinal_type,
85  class Node = typename SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal>::node_type>
86  class TekoSmoother : public SmootherPrototype<Scalar,LocalOrdinal,GlobalOrdinal,Node>
87  {
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 
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  size_t getNodeSmootherComplexity() const;
163 
165 
166  private:
167 
169  std::string type_;
170  }; // class TekoSmoother
171 
172 
187  template <class GlobalOrdinal,
188  class Node>
189  class TekoSmoother<double,int,GlobalOrdinal,Node> : public SmootherPrototype<double,int,GlobalOrdinal,Node>
190  {
191  typedef int LocalOrdinal;
192  typedef double Scalar;
194 
195 #undef MUELU_TEKOSMOOTHER_SHORT
196 #include "MueLu_UseShortNames.hpp"
197 
198  public:
199 
201 
202 
205  TekoSmoother() : type_("Teko smoother"), A_(Teuchos::null), bA_(Teuchos::null), bThyOp_(Teuchos::null), tekoParams_(Teuchos::null), inverseOp_(Teuchos::null) { };
206 
208  virtual ~TekoSmoother() { }
210 
212 
214  RCP<ParameterList> validParamList = rcp(new ParameterList());
215 
216  validParamList->set< RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A");
217  validParamList->set< std::string > ("Inverse Type", "", "Name of parameter list within 'Teko parameters' containing the Teko smoother parameters.");
218 
219  return validParamList;
220  }
221 
222  void DeclareInput(Level &currentLevel) const {
223  this->Input(currentLevel, "A");
224  }
225 
226  void SetTekoParameters(RCP<ParameterList> tekoParams) { tekoParams_ = tekoParams; };
228 
230 
231 
234  void Setup(Level &currentLevel) {
235  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
236 
237  FactoryMonitor m(*this, "Setup TekoSmoother", currentLevel);
238  if (this->IsSetup() == true)
239  this->GetOStream(Warnings0) << "MueLu::TekoSmoother::Setup(): Setup() has already been called";
240 
241  // extract blocked operator A from current level
242  A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A"); // A needed for extracting map extractors
243  bA_ = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
245  "MueLu::TekoSmoother::Build: input matrix A is not of type BlockedCrsMatrix.");
246 
247  bThyOp_ = bA_->getThyraOperator();
249  "MueLu::TekoSmoother::Build: Could not extract thyra operator from BlockedCrsMatrix.");
250 
251  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > thyOp = Teuchos::rcp_dynamic_cast<const Thyra::LinearOpBase<Scalar> >(bThyOp_);
253  "MueLu::TekoSmoother::Build: Downcast of Thyra::BlockedLinearOpBase to Teko::LinearOp failed.");
254 
255  // parameter list contains TekoSmoother parameters but does not handle the Teko parameters itself!
257  std::string smootherType = pL.get<std::string>("Inverse Type");
259  "MueLu::TekoSmoother::Build: You must provide a 'Smoother Type' name that is defined in the 'Teko parameters' sublist.");
260  type_ = smootherType;
261 
262  TEUCHOS_TEST_FOR_EXCEPTION(tekoParams_.is_null(), Exceptions::BadCast,
263  "MueLu::TekoSmoother::Build: No Teko parameters have been set.");
264 
265  Teuchos::RCP<Teko::InverseLibrary> invLib = Teko::InverseLibrary::buildFromParameterList(*tekoParams_);
266  Teuchos::RCP<Teko::InverseFactory> inverse = invLib->getInverseFactory(smootherType);
267 
268  inverseOp_ = Teko::buildInverse(*inverse, thyOp);
269  TEUCHOS_TEST_FOR_EXCEPTION(inverseOp_.is_null(), Exceptions::BadCast,
270  "MueLu::TekoSmoother::Build: Failed to build Teko inverse operator. Probably a problem with the Teko parameters.");
271 
272  this->IsSetup(true);
273  }
274 
281  void Apply(MultiVector& X, const MultiVector& B, bool /* InitialGuessIsZero */ = false) const {
283  "MueLu::TekoSmoother::Apply(): Setup() has not been called");
284 
285  Teuchos::RCP<const Teuchos::Comm<int> > comm = X.getMap()->getComm();
286 
287  Teuchos::RCP<const MapExtractor> rgMapExtractor = bA_->getRangeMapExtractor();
288  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgMapExtractor));
289 
290  // copy initial solution vector X to Ptr<Thyra::MultiVectorBase> YY
291 
292  // create a Thyra RHS vector
293  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > thyB = Thyra::createMembers(Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar> >(bThyOp_->productRange()),Teuchos::as<int>(B.getNumVectors()));
295  Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar> >(thyB);
297  "MueLu::TekoSmoother::Apply: Failed to cast range space to product range space.");
298 
299  // copy RHS vector B to Thyra::MultiVectorBase thyProdB
300  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::updateThyra(Teuchos::rcpFromRef(B), rgMapExtractor, thyProdB);
301 
302  // create a Thyra SOL vector
303  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > thyX = Thyra::createMembers(Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar> >(bThyOp_->productDomain()),Teuchos::as<int>(X.getNumVectors()));
305  Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar> >(thyX);
307  "MueLu::TekoSmoother::Apply: Failed to cast domain space to product domain space.");
308 
309  // copy RHS vector X to Thyra::MultiVectorBase thyProdX
310  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::updateThyra(Teuchos::rcpFromRef(X), rgMapExtractor, thyProdX);
311 
312  inverseOp_->apply(
313  Thyra::NOTRANS,
314  *thyB, //const MultiVectorBase<Scalar> &X,
315  thyX.ptr(), //const Ptr<MultiVectorBase<Scalar> > &Y,
316  1.0,
317  0.0);
318 
319  // copy back content of Ptr<Thyra::MultiVectorBase> thyX into X
321  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyX, comm);
322 
324  }
326 
328 
330 
331 
333  std::string description() const {
334  std::ostringstream out;
336  out << "{type = " << type_ << "}";
337  return out.str();
338  }
339 
341  //using MueLu::Describable::describe; // overloading, not hiding
342  void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel = Default) const {
344 
345  if (verbLevel & Parameters0)
346  out0 << "Prec. type: " << type_ << std::endl;
347 
348  if (verbLevel & Debug)
349  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
350  }
351 
353  size_t getNodeSmootherComplexity() const {size_t cplx=0; return cplx;}
354 
356 
357  private:
358 
360  std::string type_;
361 
363  RCP<Matrix> A_; // < ! internal blocked operator "A" generated by AFact_
366 
368  RCP<ParameterList> tekoParams_; // < ! parameter list containing Teko parameters. These parameters are not administrated by the factory and not validated.
369 
370  Teko::LinearOp inverseOp_; // < ! Teko inverse operator
371  }; // class TekoSmoother (specialization on SC=double)
372 } // namespace MueLu
373 
374 #define MUELU_TEKOSMOOTHER_SHORT
375 
376 #endif // HAVE_MUELU_TEKO
377 
378 #endif /* MUELU_TEKOSMOOTHER_DECL_HPP_ */
Important warning messages (one line)
virtual const Teuchos::ParameterList & GetParameterList() const
std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
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
RCP< const ParameterList > GetValidParameterList() const
Input.
bool is_null(const std::shared_ptr< T > &p)
RCP< SmootherPrototype > Copy() const
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual ~TekoSmoother()
Destructor.
Base class for smoother prototypes.
Print additional debugging information.
std::string description() const
Return a simple one-line description of this object.
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)
RCP< const ParameterList > GetValidParameterList() const
Input.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Ptr< T > ptr() const
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
bool IsSetup() const
Get the state of a smoother prototype.
std::string type_
smoother type
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Print class parameters.
void DeclareInput(Level &currentLevel) const
Input.
void Setup(Level &currentLevel)
Setup routine.
void Apply(MultiVector &X, const MultiVector &B, bool=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
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Teko smoother.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
virtual std::string description() const
Return a simple one-line description of this object.
std::string toString(const T &t)
int inverse(const Epetra_SerialDenseMatrix &, Epetra_SerialDenseMatrix &)
Interface to block smoothers in Teko package.
bool is_null() const