MueLu Version of the Day
Loading...
Searching...
No Matches
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
71namespace MueLu {
72
82 template <class Scalar = SmootherPrototype<>::scalar_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
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 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;
193 typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> MapExtractorClass;
194
195#undef MUELU_TEKOSMOOTHER_SHORT
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
213 RCP<const ParameterList> GetValidParameterList() const {
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_);
244 TEUCHOS_TEST_FOR_EXCEPTION(bA_.is_null(), Exceptions::BadCast,
245 "MueLu::TekoSmoother::Build: input matrix A is not of type BlockedCrsMatrix.");
246
247 bThyOp_ = bA_->getThyraOperator();
248 TEUCHOS_TEST_FOR_EXCEPTION(bThyOp_.is_null(), Exceptions::BadCast,
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_);
252 TEUCHOS_TEST_FOR_EXCEPTION(thyOp.is_null(), Exceptions::BadCast,
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!
256 const ParameterList& pL = Factory::GetParameterList();
257 std::string smootherType = pL.get<std::string>("Inverse Type");
258 TEUCHOS_TEST_FOR_EXCEPTION(smootherType.empty(), Exceptions::RuntimeError,
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 {
282 TEUCHOS_TEST_FOR_EXCEPTION(this->IsSetup() == false, Exceptions::RuntimeError,
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()));
294 Teuchos::RCP<Thyra::ProductMultiVectorBase<Scalar> > thyProdB =
295 Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar> >(thyB);
296 TEUCHOS_TEST_FOR_EXCEPTION(thyProdB.is_null(), Exceptions::BadCast,
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()));
304 Teuchos::RCP<Thyra::ProductMultiVectorBase<Scalar> > thyProdX =
305 Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar> >(thyX);
306 TEUCHOS_TEST_FOR_EXCEPTION(thyProdX.is_null(), Exceptions::BadCast,
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
320 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > XX =
321 Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyX, comm);
322
323 X.update(Teuchos::ScalarTraits<Scalar>::one(), *XX, Teuchos::ScalarTraits<Scalar>::zero());
324 }
326
327 RCP<SmootherPrototype> Copy() const { return Teuchos::rcp (new MueLu::TekoSmoother<double,int,GlobalOrdinal,Node> (*this)); }
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_
364 RCP<BlockedCrsMatrix> bA_;
365 RCP<const Thyra::BlockedLinearOpBase<Scalar> > bThyOp_;
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_ */
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
virtual std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
void Input(Level &level, const std::string &varName) const
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
virtual const Teuchos::ParameterList & GetParameterList() const =0
Base class for smoother prototypes.
bool IsSetup() const
Get the state of a smoother prototype.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
RCP< const Thyra::BlockedLinearOpBase< Scalar > > bThyOp_
std::string description() const
Return a simple one-line description of this object.
RCP< const ParameterList > GetValidParameterList() const
Input.
void Apply(MultiVector &X, const MultiVector &B, bool=false) const
Apply the Teko smoother.
Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > MapExtractorClass
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
Interface to block smoothers in Teko package.
Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > MapExtractorClass
std::string description() const
Return a simple one-line description of this object.
virtual ~TekoSmoother()
Destructor.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
RCP< const ParameterList > GetValidParameterList() const
Input.
std::string type_
smoother type
void Setup(Level &currentLevel)
Setup routine.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Teko smoother.
void SetTekoParameters(RCP< ParameterList > tekoParams)
RCP< SmootherPrototype > Copy() const
void DeclareInput(Level &currentLevel) const
Input.
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
int inverse(const Epetra_SerialDenseMatrix &, Epetra_SerialDenseMatrix &)
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.