MueLu  Version of the Day
MueLu_IfpackSmoother.cpp
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 #include "MueLu_ConfigDefs.hpp"
47 
48 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_IFPACK)
49 #include <Ifpack.h>
50 #include <Ifpack_Chebyshev.h>
51 #include "Xpetra_MultiVectorFactory.hpp"
52 
53 #include "MueLu_IfpackSmoother.hpp"
54 
55 #include "MueLu_Level.hpp"
56 #include "MueLu_Utilities.hpp"
57 #include "MueLu_Monitor.hpp"
58 
59 namespace MueLu {
60 
61  template <class Node>
62  IfpackSmoother<Node>::IfpackSmoother(std::string const & type, Teuchos::ParameterList const & paramList, LO const &overlap)
63  : type_(type), overlap_(overlap)
64  {
65  this->declareConstructionOutcome(false, "");
66  SetParameterList(paramList);
67  }
68 
69  template <class Node>
70  void IfpackSmoother<Node>::SetParameterList(const Teuchos::ParameterList& paramList) {
71  Factory::SetParameterList(paramList);
72 
74  // It might be invalid to change parameters after the setup, but it depends entirely on Ifpack implementation.
75  // TODO: I don't know if Ifpack returns an error code or exception or ignore parameters modification in this case...
76  prec_->SetParameters(const_cast<ParameterList&>(this->GetParameterList()));
77  }
78  }
79 
80  template <class Node>
81  void IfpackSmoother<Node>::SetPrecParameters(const Teuchos::ParameterList& list) const {
82  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
83  paramList.setParameters(list);
84 
85  RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
86 
87  prec_->SetParameters(*precList);
88 
89  // We would like to have the following line here:
90  // paramList.setParameters(*precList);
91  // For instance, if Ifpack sets somem parameters internally, we would like to have
92  // them listed when we call this->GetParameterList()
93  // But because of the way Ifpack handles the list, we cannot do that.
94  // The bad scenario goes like this:
95  // * SmootherFactory calls Setup
96  // * Setup calls SetPrecParameters
97  // * We call prec_->SetParameters(*precList)
98  // This actually updates the internal parameter list with default prec_ parameters
99  // This means that we get a parameter ("chebyshev: max eigenvalue", -1) in the list
100  // * Setup calls prec_->Compute()
101  // Here we may compute the max eigenvalue, but we get no indication of this. If we
102  // do compute it, our parameter list becomes outdated
103  // * SmootherFactory calls Apply
104  // * Apply constructs a list with a list with an entry "chebyshev: zero starting solution"
105  // * We call prec_->SetParameters(*precList)
106  // The last call is the problem. At this point, we have a list with an outdated entry
107  // "chebyshev: max eigenvalue", but prec_ uses this entry and replaces the computed max
108  // eigenvalue with the one from the list, resulting in -1.0 eigenvalue.
109  //
110  // Ifpack2 does not have this problem, as it does not populate the list with new entries
111  }
112 
113  template <class Node>
114  void IfpackSmoother<Node>::DeclareInput(Level &currentLevel) const {
115  this->Input(currentLevel, "A");
116 
117  if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
118  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
119  type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
120  type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
121  type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
122  type_ == "LINESMOOTHING_BLOCKRELAXATION") {
123  this->Input(currentLevel, "CoarseNumZLayers"); // necessary for fallback criterion
124  this->Input(currentLevel, "LineDetection_VertLineIds"); // necessary to feed block smoother
125  } // if (type_ == "LINESMOOTHING_BANDEDRELAXATION")
126  }
127 
128  template <class Node>
129  void IfpackSmoother<Node>::Setup(Level &currentLevel) {
130  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
131  if (SmootherPrototype::IsSetup() == true)
132  this->GetOStream(Warnings0) << "MueLu::IfpackSmoother::Setup(): Setup() has already been called" << std::endl;
133 
134  A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");
135 
136  double lambdaMax = -1.0;
137  if (type_ == "Chebyshev") {
138  std::string maxEigString = "chebyshev: max eigenvalue";
139  std::string eigRatioString = "chebyshev: ratio eigenvalue";
140 
141  try {
142  lambdaMax = Teuchos::getValue<Scalar>(this->GetParameter(maxEigString));
143  this->GetOStream(Statistics1) << maxEigString << " (cached with smoother parameter list) = " << lambdaMax << std::endl;
144 
145  } catch (Teuchos::Exceptions::InvalidParameterName&) {
146  lambdaMax = A_->GetMaxEigenvalueEstimate();
147 
148  if (lambdaMax != -1.0) {
149  this->GetOStream(Statistics1) << maxEigString << " (cached with matrix) = " << lambdaMax << std::endl;
150  this->SetParameter(maxEigString, ParameterEntry(lambdaMax));
151  }
152  }
153 
154  // Calculate the eigenvalue ratio
155  const Scalar defaultEigRatio = 20;
156 
157  Scalar ratio = defaultEigRatio;
158  try {
159  ratio = Teuchos::getValue<Scalar>(this->GetParameter(eigRatioString));
160 
161  } catch (Teuchos::Exceptions::InvalidParameterName&) {
162  this->SetParameter(eigRatioString, ParameterEntry(ratio));
163  }
164 
165  if (currentLevel.GetLevelID()) {
166  // Update ratio to be
167  // ratio = max(number of fine DOFs / number of coarse DOFs, defaultValue)
168  //
169  // NOTE: We don't need to request previous level matrix as we know for sure it was constructed
170  RCP<const Matrix> fineA = currentLevel.GetPreviousLevel()->Get<RCP<Matrix> >("A");
171  size_t nRowsFine = fineA->getGlobalNumRows();
172  size_t nRowsCoarse = A_->getGlobalNumRows();
173 
174  ratio = std::max(ratio, as<Scalar>(nRowsFine)/nRowsCoarse);
175 
176  this->GetOStream(Statistics1) << eigRatioString << " (computed) = " << ratio << std::endl;
177  this->SetParameter(eigRatioString, ParameterEntry(ratio));
178  }
179  } // if (type_ == "Chebyshev")
180 
181  if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
182  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
183  type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
184  type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
185  type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
186  type_ == "LINESMOOTHING_BLOCKRELAXATION" ) {
187  ParameterList& myparamList = const_cast<ParameterList&>(this->GetParameterList());
188 
189  LO CoarseNumZLayers = currentLevel.Get<LO>("CoarseNumZLayers",Factory::GetFactory("CoarseNumZLayers").get());
190  if (CoarseNumZLayers > 0) {
191  Teuchos::ArrayRCP<LO> TVertLineIdSmoo = currentLevel.Get< Teuchos::ArrayRCP<LO> >("LineDetection_VertLineIds", Factory::GetFactory("LineDetection_VertLineIds").get());
192 
193  // determine number of local parts
194  LO maxPart = 0;
195  for(size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
196  if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
197  }
198 
199  size_t numLocalRows = A_->getNodeNumRows();
200  TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % TVertLineIdSmoo.size() != 0, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
201 
202  if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
203  myparamList.set("partitioner: type","user");
204  myparamList.set("partitioner: map",&(TVertLineIdSmoo[0]));
205  myparamList.set("partitioner: local parts",maxPart+1);
206  } else {
207  // we assume a constant number of DOFs per node
208  size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
209 
210  // Create a new Teuchos::ArrayRCP<LO> of size numLocalRows and fill it with the corresponding information
211  Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
212  for (size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
213  for (size_t dof = 0; dof < numDofsPerNode; dof++)
214  partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
215  myparamList.set("partitioner: type","user");
216  myparamList.set("partitioner: map",&(partitionerMap[0]));
217  myparamList.set("partitioner: local parts",maxPart + 1);
218  }
219 
220  if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
221  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
222  type_ == "LINESMOOTHING_BANDEDRELAXATION")
223  type_ = "block relaxation";
224  else
225  type_ = "block relaxation";
226  } else {
227  // line detection failed -> fallback to point-wise relaxation
228  this->GetOStream(Runtime0) << "Line detection failed: fall back to point-wise relaxation" << std::endl;
229  myparamList.remove("partitioner: type",false);
230  myparamList.remove("partitioner: map", false);
231  myparamList.remove("partitioner: local parts",false);
232  type_ = "point relaxation stand-alone";
233  }
234 
235  } // if (type_ == "LINESMOOTHING_BANDEDRELAXATION")
236 
237  // If we're using a linear partitioner and haven't set the # local parts, set it to match the operator's block size
238  ParameterList precList = this->GetParameterList();
239  if(precList.isParameter("partitioner: type") && precList.get<std::string>("partitioner: type") == "linear" &&
240  !precList.isParameter("partitioner: local parts")) {
241  precList.set("partitioner: local parts", (int)A_->getNodeNumRows() / A_->GetFixedBlockSize());
242  }
243 
244 
245  RCP<Epetra_CrsMatrix> epA = Utilities::Op2NonConstEpetraCrs(A_);
246 
247  Ifpack factory;
248  prec_ = rcp(factory.Create(type_, &(*epA), overlap_));
249  TEUCHOS_TEST_FOR_EXCEPTION(prec_.is_null(), Exceptions::RuntimeError, "Could not create an Ifpack preconditioner with type = \"" << type_ << "\"");
250  SetPrecParameters();
251  prec_->Compute();
252 
254 
255  if (type_ == "Chebyshev" && lambdaMax == -1.0) {
256  Teuchos::RCP<Ifpack_Chebyshev> chebyPrec = rcp_dynamic_cast<Ifpack_Chebyshev>(prec_);
257  if (chebyPrec != Teuchos::null) {
258  lambdaMax = chebyPrec->GetLambdaMax();
259  A_->SetMaxEigenvalueEstimate(lambdaMax);
260  this->GetOStream(Statistics1) << "chebyshev: max eigenvalue (calculated by Ifpack)" << " = " << lambdaMax << std::endl;
261  }
262  TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == -1.0, Exceptions::RuntimeError, "MueLu::IfpackSmoother::Setup(): no maximum eigenvalue estimate");
263  }
264 
265  this->GetOStream(Statistics1) << description() << std::endl;
266  }
267 
268  template <class Node>
269  void IfpackSmoother<Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
270  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::IfpackSmoother::Apply(): Setup() has not been called");
271 
272  // Forward the InitialGuessIsZero option to Ifpack
273  Teuchos::ParameterList paramList;
274  bool supportInitialGuess = false;
275  if (type_ == "Chebyshev") {
276  paramList.set("chebyshev: zero starting solution", InitialGuessIsZero);
277  supportInitialGuess = true;
278 
279  } else if (type_ == "point relaxation stand-alone") {
280  paramList.set("relaxation: zero starting solution", InitialGuessIsZero);
281  supportInitialGuess = true;
282  }
283 
284  SetPrecParameters(paramList);
285 
286  // Apply
287  if (InitialGuessIsZero || supportInitialGuess) {
290 
291  prec_->ApplyInverse(epB, epX);
292 
293  } else {
294  RCP<MultiVector> Residual = Utilities::Residual(*A_, X, B);
295  RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
296 
298  const Epetra_MultiVector& epB = Utilities::MV2EpetraMV(*Residual);
299 
300  prec_->ApplyInverse(epB, epX);
301 
302  X.update(1.0, *Correction, 1.0);
303  }
304  }
305 
306  template <class Node>
307  RCP<MueLu::SmootherPrototype<double, int, int, Node> > IfpackSmoother<Node>::Copy() const {
308  RCP<IfpackSmoother<Node> > smoother = rcp(new IfpackSmoother<Node>(*this) );
309  smoother->SetParameterList(this->GetParameterList());
310  return Teuchos::rcp_dynamic_cast<MueLu::SmootherPrototype<double, int, int, Node> >(smoother);
311  }
312 
313  template <class Node>
314  std::string IfpackSmoother<Node>::description() const {
315  std::ostringstream out;
316  // The check "GetVerbLevel() == Test" is to avoid
317  // failures in the EasyInterface test.
318  if (prec_ == Teuchos::null || this->GetVerbLevel() == Test) {
320  out << "{type = " << type_ << "}";
321  } else {
322  out << prec_->Label();
323  }
324  return out.str();
325  }
326 
327  template <class Node>
328  void IfpackSmoother<Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
330 
331  if (verbLevel & Parameters0)
332  out0 << "Prec. type: " << type_ << std::endl;
333 
334  if (verbLevel & Parameters1) {
335  out0 << "Parameter list: " << std::endl;
336  Teuchos::OSTab tab2(out);
337  out << this->GetParameterList();
338  out0 << "Overlap: " << overlap_ << std::endl;
339  }
340 
341  if (verbLevel & External)
342  if (prec_ != Teuchos::null) {
343  Teuchos::OSTab tab2(out);
344  out << *prec_ << std::endl;
345  }
346 
347  if (verbLevel & Debug) {
348  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
349  << "-" << std::endl
350  << "RCP<A_>: " << A_ << std::endl
351  << "RCP<prec_>: " << prec_ << std::endl;
352  }
353  }
354 
355  template <class Node>
357  // FIXME: This is a placeholder
358  return Teuchos::OrdinalTraits<size_t>::invalid();
359  }
360 
361 
362 } // namespace MueLu
363 
364 // The IfpackSmoother is only templated on the Node, since it is an Epetra only object
365 // Therefore we do not need the full ETI instantiations as we do for the other MueLu
366 // objects which are instantiated on all template parameters.
367 #if defined(HAVE_MUELU_EPETRA)
369 #endif
370 
371 #endif
void Setup(Level &currentLevel)
Set up the smoother.
Important warning messages (one line)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Print external lib objects.
bool IsSetup() const
Get the state of a smoother prototype.
static RCP< const Epetra_MultiVector > MV2EpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
Timer to be used in factories. Similar to Monitor but with additional timers.
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
Print more statistics.
RCP< SmootherPrototype > Copy() const
Print additional debugging information.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
Ifpack_Preconditioner * Create(const std::string PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
IfpackSmoother(std::string const &type, Teuchos::ParameterList const &paramList=Teuchos::ParameterList(), LO const &overlap=0)
Constructor.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
void declareConstructionOutcome(bool fail, std::string msg)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
virtual double GetLambdaMax()
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
void DeclareInput(Level &currentLevel) const
Input.
void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
Print class parameters.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
std::string description() const
Return a simple one-line description of this object.
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
virtual std::string description() const
Return a simple one-line description of this object.
Class that encapsulates Ifpack smoothers.