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