MueLu  Version of the Day
MueLu_Ifpack2Smoother_def.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 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_IFPACK2SMOOTHER_DEF_HPP
47 #define MUELU_IFPACK2SMOOTHER_DEF_HPP
48 
49 #include "MueLu_ConfigDefs.hpp"
50 
51 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2)
52 
53 #include <Teuchos_ParameterList.hpp>
54 
55 #include <Tpetra_RowMatrix.hpp>
56 
57 #include <Ifpack2_Chebyshev.hpp>
58 #include <Ifpack2_Factory.hpp>
59 #include <Ifpack2_Parameters.hpp>
60 
61 #include <Xpetra_BlockedCrsMatrix.hpp>
62 #include <Xpetra_CrsMatrix.hpp>
63 #include <Xpetra_CrsMatrixWrap.hpp>
64 #include <Xpetra_Matrix.hpp>
65 #include <Xpetra_MultiVectorFactory.hpp>
66 
68 #include "MueLu_Level.hpp"
70 #include "MueLu_Utilities.hpp"
71 #include "MueLu_Monitor.hpp"
72 
73 // #define IFPACK2_HAS_PROPER_REUSE
74 
75 namespace MueLu {
76 
77  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
78  Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Ifpack2Smoother(const std::string& type, const Teuchos::ParameterList& paramList, const LO& overlap)
79  : type_(type), overlap_(overlap)
80  {
81  SetParameterList(paramList);
82  }
83 
84  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
86  Factory::SetParameterList(paramList);
87 
89  // It might be invalid to change parameters after the setup, but it depends entirely on Ifpack implementation.
90  // TODO: I don't know if Ifpack returns an error code or exception or ignore parameters modification in this case...
91  SetPrecParameters();
92  }
93  }
94 
95  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
97  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
98  paramList.setParameters(list);
99 
100  RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
101 
102  prec_->setParameters(*precList);
103 
104  paramList.setParameters(*precList); // what about that??
105  }
106 
107  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
109  this->Input(currentLevel, "A");
110 
111  if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
112  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
113  type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
114  type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
115  type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
116  type_ == "LINESMOOTHING_BLOCKRELAXATION") {
117  this->Input(currentLevel, "CoarseNumZLayers"); // necessary for fallback criterion
118  this->Input(currentLevel, "LineDetection_VertLineIds"); // necessary to feed block smoother
119  }
120  }
121 
122  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
124  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
125 
126  A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");
127 
128  if (type_ == "SCHWARZ")
129  SetupSchwarz(currentLevel);
130 
131  else if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
132  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
133  type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
134  type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
135  type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
136  type_ == "LINESMOOTHING_BLOCKRELAXATION")
137  SetupLineSmoothing(currentLevel);
138 
139  else if (type_ == "CHEBYSHEV")
140  SetupChebyshev(currentLevel);
141 
142  else
143  SetupGeneric(currentLevel);
144 
146 
147  this->GetOStream(Statistics1) << description() << std::endl;
148  }
149 
150  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
152  typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
153 
154  bool reusePreconditioner = false;
155  if (this->IsSetup() == true) {
156  // Reuse the constructed preconditioner
157  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupSchwarz(): Setup() has already been called, assuming reuse" << std::endl;
158 
159  bool isTRowMatrix = true;
160  RCP<const tRowMatrix> tA;
161  try {
163  } catch (Exceptions::BadCast) {
164  isTRowMatrix = false;
165  }
166 
167  RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
168  if (!prec.is_null() && isTRowMatrix) {
169 #ifdef IFPACK2_HAS_PROPER_REUSE
170  prec->resetMatrix(tA);
171  reusePreconditioner = true;
172 #else
173  this->GetOStream(Errors) << "Ifpack2 does not have proper reuse yet." << std::endl;
174 #endif
175 
176  } else {
177  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available "
178  "(either failed cast to CanChangeMatrix, or to Tpetra Row Matrix), reverting to full construction" << std::endl;
179  }
180  }
181 
182  if (!reusePreconditioner) {
183  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
184 
185  bool isBlockedMatrix = false;
186  RCP<Matrix> merged2Mat;
187 
188  std::string sublistName = "subdomain solver parameters";
189  if (paramList.isSublist(sublistName)) {
190  // If we are doing "user" partitioning, we assume that what the user
191  // really wants to do is make tiny little subdomains with one row
192  // assigned to each subdomain. The rows used for these little
193  // subdomains correspond to those in the 2nd block row. Then,
194  // if we overlap these mini-subdomains, we will do something that
195  // looks like Vanka (grabbing all velocities associated with each
196  // each pressure unknown). In addition, we put all Dirichlet points
197  // as a little mini-domain.
198  ParameterList& subList = paramList.sublist(sublistName);
199 
200  std::string partName = "partitioner: type";
201  if (subList.isParameter(partName) && subList.get<std::string>(partName) == "user") {
202  isBlockedMatrix = true;
203 
204  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
205  TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
206  "Matrix A must be of type BlockedCrsMatrix.");
207 
208  size_t numVels = bA->getMatrix(0,0)->getNodeNumRows();
209  size_t numPres = bA->getMatrix(1,0)->getNodeNumRows();
210  size_t numRows = A_->getNodeNumRows();
211 
212  ArrayRCP<LocalOrdinal> blockSeeds(numRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
213 
214  size_t numBlocks = 0;
215  for (size_t rowOfB = numVels; rowOfB < numVels+numPres; ++rowOfB)
216  blockSeeds[rowOfB] = numBlocks++;
217 
218  RCP<BlockedCrsMatrix> bA2 = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
219  TEUCHOS_TEST_FOR_EXCEPTION(bA2.is_null(), Exceptions::BadCast,
220  "Matrix A must be of type BlockedCrsMatrix.");
221 
222  merged2Mat = bA2->Merge();
223 
224  // Add Dirichlet rows to the list of seeds
225  ArrayRCP<const bool> boundaryNodes = Utilities::DetectDirichletRows(*merged2Mat, 0.0);
226  bool haveBoundary = false;
227  for (LO i = 0; i < boundaryNodes.size(); i++)
228  if (boundaryNodes[i]) {
229  // FIXME:
230  // 1. would not this [] overlap with some in the previos blockSeed loop?
231  // 2. do we need to distinguish between pressure and velocity Dirichlet b.c.
232  blockSeeds[i] = numBlocks;
233  haveBoundary = true;
234  }
235  if (haveBoundary)
236  numBlocks++;
237 
238  subList.set("partitioner: map", blockSeeds);
239  subList.set("partitioner: local parts", as<int>(numBlocks));
240 
241  } else {
242  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
243  if (!bA.is_null()) {
244  isBlockedMatrix = true;
245  merged2Mat = bA->Merge();
246  }
247  }
248  }
249 
250  RCP<const tRowMatrix> tA;
251  if (isBlockedMatrix == true) tA = Utilities::Op2NonConstTpetraRow(merged2Mat);
252  else tA = Utilities::Op2NonConstTpetraRow(A_);
253 
254  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
255  SetPrecParameters();
256 
257  prec_->initialize();
258  }
259 
260  prec_->compute();
261  }
262 
263  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
265  if (this->IsSetup() == true) {
266  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupLineSmoothing(): Setup() has already been called" << std::endl;
267  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupLineSmoothing(): reuse of this type is not available, reverting to full construction" << std::endl;
268  }
269 
270  ParameterList& myparamList = const_cast<ParameterList&>(this->GetParameterList());
271 
272  LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,"CoarseNumZLayers");
273  if (CoarseNumZLayers > 0) {
274  Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get< Teuchos::ArrayRCP<LO> >(currentLevel, "LineDetection_VertLineIds");
275 
276  // determine number of local parts
277  LO maxPart = 0;
278  for(size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
279  if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
280  }
281 
282  size_t numLocalRows = A_->getNodeNumRows();
283  TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % TVertLineIdSmoo.size() != 0, Exceptions::RuntimeError,
284  "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
285 
286  if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
287  myparamList.set("partitioner: type","user");
288  myparamList.set("partitioner: map",TVertLineIdSmoo);
289  myparamList.set("partitioner: local parts",maxPart+1);
290  } else {
291  // we assume a constant number of DOFs per node
292  size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
293 
294  // Create a new Teuchos::ArrayRCP<LO> of size numLocalRows and fill it with the corresponding information
295  Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
296  for (size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
297  for (size_t dof = 0; dof < numDofsPerNode; dof++)
298  partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
299  myparamList.set("partitioner: type","user");
300  myparamList.set("partitioner: map",partitionerMap);
301  myparamList.set("partitioner: local parts",maxPart + 1);
302  }
303 
304  if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
305  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
306  type_ == "LINESMOOTHING_BANDEDRELAXATION")
307  type_ = "BANDEDRELAXATION";
308  else
309  type_ = "BLOCKRELAXATION";
310  } else {
311  // line detection failed -> fallback to point-wise relaxation
312  this->GetOStream(Runtime0) << "Line detection failed: fall back to point-wise relaxation" << std::endl;
313  myparamList.remove("partitioner: type",false);
314  myparamList.remove("partitioner: map", false);
315  myparamList.remove("partitioner: local parts",false);
316  type_ = "RELAXATION";
317  }
318 
319  RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tA = Utilities::Op2NonConstTpetraRow(A_);
320 
321  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
322  SetPrecParameters();
323  prec_->initialize();
324  prec_->compute();
325  }
326 
327  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
329  if (this->IsSetup() == true) {
330  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupChebyshev(): SetupChebyshev() has already been called" << std::endl;
331  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupChebyshev(): reuse of this type is not available, reverting to full construction" << std::endl;
332  }
333 
334  typedef Teuchos::ScalarTraits<SC> STS;
335  SC negone = -STS::one();
336 
337  SC lambdaMax = negone;
338  {
339  std::string maxEigString = "chebyshev: max eigenvalue";
340  std::string eigRatioString = "chebyshev: ratio eigenvalue";
341 
342  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
343 
344  // Get/calculate the maximum eigenvalue
345  if (paramList.isParameter(maxEigString)) {
346  if (paramList.isType<double>(maxEigString))
347  lambdaMax = paramList.get<double>(maxEigString);
348  else
349  lambdaMax = paramList.get<SC>(maxEigString);
350  this->GetOStream(Statistics1) << maxEigString << " (cached with smoother parameter list) = " << lambdaMax << std::endl;
351 
352  } else {
353  lambdaMax = A_->GetMaxEigenvalueEstimate();
354  if (lambdaMax != negone) {
355  this->GetOStream(Statistics1) << maxEigString << " (cached with matrix) = " << lambdaMax << std::endl;
356  paramList.set(maxEigString, lambdaMax);
357  }
358  }
359 
360  // Calculate the eigenvalue ratio
361  const SC defaultEigRatio = 20;
362 
363  SC ratio = defaultEigRatio;
364  if (paramList.isParameter(eigRatioString)) {
365  if (paramList.isType<double>(eigRatioString))
366  ratio = paramList.get<double>(eigRatioString);
367  else
368  ratio = paramList.get<SC>(eigRatioString);
369  }
370  if (currentLevel.GetLevelID()) {
371  // Update ratio to be
372  // ratio = max(number of fine DOFs / number of coarse DOFs, defaultValue)
373  //
374  // NOTE: We don't need to request previous level matrix as we know for sure it was constructed
375  RCP<const Matrix> fineA = currentLevel.GetPreviousLevel()->Get<RCP<Matrix> >("A");
376  size_t nRowsFine = fineA->getGlobalNumRows();
377  size_t nRowsCoarse = A_->getGlobalNumRows();
378 
379  SC levelRatio = as<SC>(as<float>(nRowsFine)/nRowsCoarse);
380  if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
381  ratio = levelRatio;
382  }
383 
384  this->GetOStream(Statistics1) << eigRatioString << " (computed) = " << ratio << std::endl;
385  paramList.set(eigRatioString, ratio);
386  }
387 
388  RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tA = Utilities::Op2NonConstTpetraRow(A_);
389 
390  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
391  SetPrecParameters();
392  prec_->initialize();
393  prec_->compute();
394 
395  if (lambdaMax == negone) {
396  typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
397 
398  Teuchos::RCP<Ifpack2::Chebyshev<MatrixType> > chebyPrec = rcp_dynamic_cast<Ifpack2::Chebyshev<MatrixType> >(prec_);
399  if (chebyPrec != Teuchos::null) {
400  lambdaMax = chebyPrec->getLambdaMaxForApply();
401  A_->SetMaxEigenvalueEstimate(lambdaMax);
402  this->GetOStream(Statistics1) << "chebyshev: max eigenvalue (calculated by Ifpack2)" << " = " << lambdaMax << std::endl;
403  }
404  TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == negone, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Setup(): no maximum eigenvalue estimate");
405  }
406  }
407 
408  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
410  typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
411 
412  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
413  if (!bA.is_null())
414  A_ = bA->Merge();
415 
416  RCP<const tRowMatrix> tA = Utilities::Op2NonConstTpetraRow(A_);
417 
418  bool reusePreconditioner = false;
419  if (this->IsSetup() == true) {
420  // Reuse the constructed preconditioner
421  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
422 
423  RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
424  if (!prec.is_null()) {
425 #ifdef IFPACK2_HAS_PROPER_REUSE
426  prec->resetMatrix(tA);
427  reusePreconditioner = true;
428 #else
429  this->GetOStream(Errors) << "Ifpack2 does not have proper reuse yet." << std::endl;
430 #endif
431 
432  } else {
433  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available (failed cast to CanChangeMatrix), "
434  "reverting to full construction" << std::endl;
435  }
436  }
437 
438  if (!reusePreconditioner) {
439  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
440  SetPrecParameters();
441  prec_->initialize();
442  }
443 
444  prec_->compute();
445  }
446 
447  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
448  void Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
449  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Apply(): Setup() has not been called");
450 
451  // Forward the InitialGuessIsZero option to Ifpack2
452  // TODO: It might be nice to switch back the internal
453  // "zero starting solution" option of the ifpack2 object prec_ to his
454  // initial value at the end but there is no way right now to get
455  // the current value of the "zero starting solution" in ifpack2.
456  // It's not really an issue, as prec_ can only be used by this method.
457  // TODO: When https://software.sandia.gov/bugzilla/show_bug.cgi?id=5283#c2 is done
458  // we should remove the if/else/elseif and just test if this
459  // option is supported by current ifpack2 preconditioner
460  Teuchos::ParameterList paramList;
461  bool supportInitialGuess = false;
462  if (type_ == "CHEBYSHEV") {
463  paramList.set("chebyshev: zero starting solution", InitialGuessIsZero);
464  SetPrecParameters(paramList);
465  supportInitialGuess = true;
466 
467  } else if (type_ == "RELAXATION") {
468  paramList.set("relaxation: zero starting solution", InitialGuessIsZero);
469  SetPrecParameters(paramList);
470  supportInitialGuess = true;
471 
472  } else if (type_ == "KRYLOV") {
473  paramList.set("krylov: zero starting solution", InitialGuessIsZero);
474  SetPrecParameters(paramList);
475  supportInitialGuess = true;
476 
477  } else if (type_ == "SCHWARZ") {
478  paramList.set("schwarz: zero starting solution", InitialGuessIsZero);
479  //Because additive Schwarz has "delta" semantics, it's sufficient to
480  //toggle only the zero initial guess flag, and not pass in already
481  //set parameters. If we call SetPrecParameters, the subdomain solver
482  //will be destroyed.
483  prec_->setParameters(paramList);
484  supportInitialGuess = true;
485  }
486 
487  //TODO JJH 30Apr2014 Calling SetPrecParameters(paramList) when the smoother
488  //is Ifpack2::AdditiveSchwarz::setParameterList() will destroy the subdomain
489  //(aka inner) solver. This behavior is documented but a departure from what
490  //it previously did, and what other Ifpack2 solvers currently do. So I have
491  //moved SetPrecParameters(paramList) into the if-else block above.
492 
493  // Apply
494  if (InitialGuessIsZero || supportInitialGuess) {
497  prec_->apply(tpB, tpX);
498  } else {
499  typedef Teuchos::ScalarTraits<Scalar> TST;
500  RCP<MultiVector> Residual = Utilities::Residual(*A_, X, B);
501  RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
502 
505 
506  prec_->apply(tpB, tpX);
507 
508  X.update(TST::one(), *Correction, TST::one());
509  }
510  }
511 
512  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
513  RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Copy() const {
514  RCP<Ifpack2Smoother> smoother = rcp(new Ifpack2Smoother(*this) );
515  smoother->SetParameterList(this->GetParameterList());
516  return smoother;
517  }
518 
519  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
521  std::ostringstream out;
523  out << prec_->description();
524  } else {
526  out << "{type = " << type_ << "}";
527  }
528  return out.str();
529  }
530 
531  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
532  void Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
534 
535  if (verbLevel & Parameters0)
536  out0 << "Prec. type: " << type_ << std::endl;
537 
538  if (verbLevel & Parameters1) {
539  out0 << "Parameter list: " << std::endl;
540  Teuchos::OSTab tab2(out);
541  out << this->GetParameterList();
542  out0 << "Overlap: " << overlap_ << std::endl;
543  }
544 
545  if (verbLevel & External)
546  if (prec_ != Teuchos::null) {
547  Teuchos::OSTab tab2(out);
548  out << *prec_ << std::endl;
549  }
550 
551  if (verbLevel & Debug) {
552  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
553  << "-" << std::endl
554  << "RCP<prec_>: " << prec_ << std::endl;
555  }
556  }
557 
558 } // namespace MueLu
559 
560 #endif // HAVE_MUELU_TPETRA && HAVE_MUELU_IFPACK2
561 #endif // MUELU_IFPACK2SMOOTHER_DEF_HPP
Important warning messages (one line)
void SetupGeneric(Level &currentLevel)
RCP< SmootherPrototype > Copy() const
Exception indicating invalid cast attempted.
RCP< Level > & GetPreviousLevel()
Previous level.
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
std::string description() const
Return a simple one-line description of this object.
Print external lib objects.
bool IsSetup() const
Get the state of a smoother prototype.
friend class Ifpack2Smoother
Constructor.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
Print additional debugging information.
One-liner description of what is happening.
void SetParameterList(const Teuchos::ParameterList &paramList)
Namespace for MueLu classes and methods.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
void SetupSchwarz(Level &currentLevel)
void SetupLineSmoothing(Level &currentLevel)
MatrixType::scalar_type getLambdaMaxForApply() const
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.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Setup(Level &currentLevel)
Set up the smoother.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
static Teuchos::RCP< Preconditioner< typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type > > create(const std::string &precType, const Teuchos::RCP< const MatrixType > &matrix)
void DeclareInput(Level &currentLevel) const
Input.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Print class parameters.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.))
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
void SetupChebyshev(Level &currentLevel)
Description of what is happening (more verbose)
Class that encapsulates Ifpack2 smoothers.
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
virtual std::string description() const
Return a simple one-line description of this object.