MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
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>
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>
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 
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 
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 = currentLevel.Get<LO>("CoarseNumZLayers",Factory::GetFactory("CoarseNumZLayers").get());
189  if (CoarseNumZLayers > 0) {
190  Teuchos::ArrayRCP<LO> TVertLineIdSmoo = currentLevel.Get< Teuchos::ArrayRCP<LO> >("LineDetection_VertLineIds", Factory::GetFactory("LineDetection_VertLineIds").get());
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
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 
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) {
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 
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>
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  template <class Node>
348  // FIXME: This is a placeholder
350  }
351 
352 
353 } // namespace MueLu
354 
355 // The IfpackSmoother is only templated on the Node, since it is an Epetra only object
356 // Therefore we do not need the full ETI instantiations as we do for the other MueLu
357 // objects which are instantiated on all template parameters.
358 #if defined(HAVE_MUELU_EPETRA)
360 #endif
361 
362 #endif
void Setup(Level &currentLevel)
Set up the smoother.
Important warning messages (one line)
RCP< SmootherPrototype > Copy() const
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Print external lib objects.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Print more statistics.
Print additional debugging information.
LocalOrdinal LO
One-liner description of what is happening.
size_type size() const
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.
bool remove(std::string const &name, bool throwIfNotExists=true)
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.
void DeclareInput(Level &currentLevel) const
Input.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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
bool IsSetup() const
Get the state of a smoother prototype.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
virtual double GetLambdaMax()
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
ParameterList & setParameters(const ParameterList &source)
void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
Print class parameters.
std::string description() const
Return a simple one-line description of this object.
virtual const RCP< const FactoryBase > GetFactory(const std::string &varName) const =0
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.
std::string toString(const T &t)