MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_Ifpack2Smoother_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_IFPACK2SMOOTHER_DEF_HPP
11 #define MUELU_IFPACK2SMOOTHER_DEF_HPP
12 
13 #include "MueLu_ConfigDefs.hpp"
15 
16 #if defined(HAVE_MUELU_IFPACK2)
17 
19 
20 #include <Tpetra_RowMatrix.hpp>
21 
22 #include <Ifpack2_Chebyshev.hpp>
23 #include <Ifpack2_Hiptmair.hpp>
24 #include <Ifpack2_RILUK.hpp>
25 #include <Ifpack2_Relaxation.hpp>
26 #include <Ifpack2_ILUT.hpp>
27 #include <Ifpack2_BlockRelaxation.hpp>
28 #include <Ifpack2_Factory.hpp>
29 #include <Ifpack2_Parameters.hpp>
30 
32 #include <Xpetra_CrsMatrix.hpp>
33 #include <Xpetra_CrsMatrixWrap.hpp>
34 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
35 #include <Xpetra_Matrix.hpp>
36 #include <Xpetra_MatrixMatrix.hpp>
37 #include <Xpetra_MultiVectorFactory.hpp>
38 #include <Xpetra_TpetraMultiVector.hpp>
39 
40 #include <Tpetra_BlockCrsMatrix_Helpers.hpp>
41 
43 #include "MueLu_Level.hpp"
44 #include "MueLu_Utilities.hpp"
45 #include "MueLu_Monitor.hpp"
46 #include "MueLu_Aggregates.hpp"
47 
48 #ifdef HAVE_MUELU_INTREPID2
51 #include "Intrepid2_Basis.hpp"
52 #include "Kokkos_DynRankView.hpp"
53 #endif
54 
55 namespace MueLu {
56 
57 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
59  : type_(type)
60  , overlap_(overlap) {
61  typedef Tpetra::RowMatrix<SC, LO, GO, NO> tRowMatrix;
62  bool isSupported = Ifpack2::Factory::isSupported<tRowMatrix>(type_) || (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
63  type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
64  type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
65  type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
66  type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
67  type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION" ||
68  type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
69  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
70  type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
71  type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
72  type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
73  type_ == "LINESMOOTHING_BLOCKRELAXATION" ||
74  type_ == "TOPOLOGICAL" ||
75  type_ == "AGGREGATE");
76  this->declareConstructionOutcome(!isSupported, "Ifpack2 does not provide the smoother '" + type_ + "'.");
77  if (isSupported)
78  SetParameterList(paramList);
79 }
80 
81 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
83  Factory::SetParameterList(paramList);
84 
86  // It might be invalid to change parameters after the setup, but it depends entirely on Ifpack implementation.
87  // TODO: I don't know if Ifpack returns an error code or exception or ignore parameters modification in this case...
88  SetPrecParameters();
89  }
90 }
91 
92 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
94  std::string prefix = this->ShortClassName() + ": SetPrecParameters";
95  RCP<TimeMonitor> tM = rcp(new TimeMonitor(*this, prefix, Timings0));
96  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
97 
98  paramList.setParameters(list);
99 
100  RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
101 
102  // Do we want an Ifpack2 apply timer?
103  precList->set("timer for apply", this->IsPrint(Timings));
104 
105  if (!precList.is_null() && precList->isParameter("partitioner: type") && precList->get<std::string>("partitioner: type") == "linear" &&
106  !precList->isParameter("partitioner: local parts")) {
107  LO matrixBlockSize = 1;
108  int lclSize = A_->getRangeMap()->getLocalNumElements();
109  RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
110  if (!matA.is_null()) {
111  lclSize = matA->getLocalNumRows();
112  matrixBlockSize = matA->GetFixedBlockSize();
113  }
114  precList->set("partitioner: local parts", lclSize / matrixBlockSize);
115  }
116 
117  prec_->setParameters(*precList);
118 
119  paramList.setParameters(*precList); // what about that??
120 }
121 
122 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
124  this->Input(currentLevel, "A");
125 
126  if (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
127  type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
128  type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
129  type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
130  type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
131  type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION" ||
132  type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
133  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
134  type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
135  type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
136  type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
137  type_ == "LINESMOOTHING_BLOCKRELAXATION") {
138  this->Input(currentLevel, "CoarseNumZLayers"); // necessary for fallback criterion
139  this->Input(currentLevel, "LineDetection_VertLineIds"); // necessary to feed block smoother
140  } else if (type_ == "BLOCK RELAXATION" ||
141  type_ == "BLOCK_RELAXATION" ||
142  type_ == "BLOCKRELAXATION" ||
143  // Banded
144  type_ == "BANDED_RELAXATION" ||
145  type_ == "BANDED RELAXATION" ||
146  type_ == "BANDEDRELAXATION" ||
147  // Tridiagonal
148  type_ == "TRIDI_RELAXATION" ||
149  type_ == "TRIDI RELAXATION" ||
150  type_ == "TRIDIRELAXATION" ||
151  type_ == "TRIDIAGONAL_RELAXATION" ||
152  type_ == "TRIDIAGONAL RELAXATION" ||
153  type_ == "TRIDIAGONALRELAXATION") {
154  // We need to check for the "partitioner type" = "line"
155  ParameterList precList = this->GetParameterList();
156  if (precList.isParameter("partitioner: type") &&
157  precList.get<std::string>("partitioner: type") == "line") {
158  this->Input(currentLevel, "Coordinates");
159  }
160  } else if (type_ == "TOPOLOGICAL") {
161  // for the topological smoother, we require an element to node map:
162  this->Input(currentLevel, "pcoarsen: element to node map");
163  } else if (type_ == "AGGREGATE") {
164  // Aggregate smoothing needs aggregates
165  this->Input(currentLevel, "Aggregates");
166  } else if (type_ == "HIPTMAIR") {
167  // Hiptmair needs D0 and NodeMatrix
168  this->Input(currentLevel, "NodeMatrix");
169  this->Input(currentLevel, "D0");
170  }
171 }
172 
173 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
175  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
176  A_ = Factory::Get<RCP<Operator>>(currentLevel, "A");
177  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
178 
179  // If the user asked us to convert the matrix into BlockCrsMatrix form, we do that now
180 
181  if (paramList.isParameter("smoother: use blockcrsmatrix storage") && paramList.get<bool>("smoother: use blockcrsmatrix storage")) {
182  int blocksize = 1;
183  RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
184  if (!matA.is_null())
185  blocksize = matA->GetFixedBlockSize();
186  if (blocksize) {
187  // NOTE: Don't think you can move this out of the if block. You can't. The test MueLu_MeshTyingBlocked_SimpleSmoother_2dof_medium_MPI_1 will fail
188 
189  RCP<CrsMatrixWrap> AcrsWrap = rcp_dynamic_cast<CrsMatrixWrap>(A_);
190  if (AcrsWrap.is_null())
191  throw std::runtime_error("Ifpack2Smoother: Cannot convert matrix A to CrsMatrixWrap object.");
192  RCP<CrsMatrix> Acrs = AcrsWrap->getCrsMatrix();
193  if (Acrs.is_null())
194  throw std::runtime_error("Ifpack2Smoother: Cannot extract CrsMatrix from matrix A.");
195  RCP<TpetraCrsMatrix> At = rcp_dynamic_cast<TpetraCrsMatrix>(Acrs);
196  if (At.is_null()) {
198  throw std::runtime_error("Ifpack2Smoother: Cannot extract CrsMatrix or BlockCrsMatrix from matrix A.");
199  this->GetOStream(Statistics0) << "Ifpack2Smoother: Using (native) BlockCrsMatrix storage with blocksize " << blocksize << std::endl;
200  paramList.remove("smoother: use blockcrsmatrix storage");
201  } else {
202  RCP<Tpetra::BlockCrsMatrix<Scalar, LO, GO, Node>> blockCrs = Tpetra::convertToBlockCrsMatrix(*At->getTpetra_CrsMatrix(), blocksize);
203  RCP<CrsMatrix> blockCrs_as_crs = rcp(new TpetraBlockCrsMatrix(blockCrs));
204  RCP<CrsMatrixWrap> blockWrap = rcp(new CrsMatrixWrap(blockCrs_as_crs));
205  A_ = blockWrap;
206  this->GetOStream(Statistics0) << "Ifpack2Smoother: Using BlockCrsMatrix storage with blocksize " << blocksize << std::endl;
207 
208  paramList.remove("smoother: use blockcrsmatrix storage");
209  }
210  }
211  }
212 
213  if (type_ == "SCHWARZ")
214  SetupSchwarz(currentLevel);
215 
216  else if (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
217  type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
218  type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
219  type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
220  type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
221  type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION" ||
222  type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
223  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
224  type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
225  type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
226  type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
227  type_ == "LINESMOOTHING_BLOCKRELAXATION")
228  SetupLineSmoothing(currentLevel);
229 
230  else if (type_ == "BLOCK_RELAXATION" ||
231  type_ == "BLOCK RELAXATION" ||
232  type_ == "BLOCKRELAXATION" ||
233  // Banded
234  type_ == "BANDED_RELAXATION" ||
235  type_ == "BANDED RELAXATION" ||
236  type_ == "BANDEDRELAXATION" ||
237  // Tridiagonal
238  type_ == "TRIDI_RELAXATION" ||
239  type_ == "TRIDI RELAXATION" ||
240  type_ == "TRIDIRELAXATION" ||
241  type_ == "TRIDIAGONAL_RELAXATION" ||
242  type_ == "TRIDIAGONAL RELAXATION" ||
243  type_ == "TRIDIAGONALRELAXATION")
244  SetupBlockRelaxation(currentLevel);
245 
246  else if (type_ == "CHEBYSHEV")
247  SetupChebyshev(currentLevel);
248 
249  else if (type_ == "TOPOLOGICAL") {
250 #ifdef HAVE_MUELU_INTREPID2
251  SetupTopological(currentLevel);
252 #else
253  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "'TOPOLOGICAL' smoother choice requires Intrepid2");
254 #endif
255  } else if (type_ == "AGGREGATE")
256  SetupAggregate(currentLevel);
257 
258  else if (type_ == "HIPTMAIR")
259  SetupHiptmair(currentLevel);
260 
261  else {
262  SetupGeneric(currentLevel);
263  }
264 
266 
267  this->GetOStream(Statistics1) << description() << std::endl;
268 }
269 
270 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
272  typedef Tpetra::RowMatrix<SC, LO, GO, NO> tRowMatrix;
273 
274  bool reusePreconditioner = false;
275  if (this->IsSetup() == true) {
276  // Reuse the constructed preconditioner
277  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupSchwarz(): Setup() has already been called, assuming reuse" << std::endl;
278 
279  bool isTRowMatrix = true;
281  try {
282  tA = Xpetra::toTpetraRowMatrix(A_);
283  } catch (Exceptions::BadCast&) {
284  isTRowMatrix = false;
285  }
286 
288  if (!prec.is_null() && isTRowMatrix) {
289  prec->setMatrix(tA);
290  reusePreconditioner = true;
291  } else {
292  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available "
293  "(either failed cast to CanChangeMatrix, or to Tpetra Row Matrix), reverting to full construction"
294  << std::endl;
295  }
296  }
297 
298  if (!reusePreconditioner) {
299  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
300 
301  bool isBlockedMatrix = false;
302  RCP<Matrix> merged2Mat;
303 
304  std::string sublistName = "subdomain solver parameters";
305  if (paramList.isSublist(sublistName)) {
306  // If we are doing "user" partitioning, we assume that what the user
307  // really wants to do is make tiny little subdomains with one row
308  // assigned to each subdomain. The rows used for these little
309  // subdomains correspond to those in the 2nd block row. Then,
310  // if we overlap these mini-subdomains, we will do something that
311  // looks like Vanka (grabbing all velocities associated with each
312  // each pressure unknown). In addition, we put all Dirichlet points
313  // as a little mini-domain.
314  ParameterList& subList = paramList.sublist(sublistName);
315 
316  std::string partName = "partitioner: type";
317  // Pretty sure no one has been using this. Unfortunately, old if
318  // statement (which checked for equality with "user") prevented
319  // anyone from employing other types of Ifpack2 user partition
320  // options. Leaving this and switching if to "vanka user" just
321  // in case some day someone might want to use this.
322  if (subList.isParameter(partName) && subList.get<std::string>(partName) == "vanka user") {
323  isBlockedMatrix = true;
324 
325  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
327  "Matrix A must be of type BlockedCrsMatrix.");
328 
329  size_t numVels = bA->getMatrix(0, 0)->getLocalNumRows();
330  size_t numPres = bA->getMatrix(1, 0)->getLocalNumRows();
331  size_t numRows = rcp_dynamic_cast<Matrix>(A_, true)->getLocalNumRows();
332 
334 
335  size_t numBlocks = 0;
336  for (size_t rowOfB = numVels; rowOfB < numVels + numPres; ++rowOfB)
337  blockSeeds[rowOfB] = numBlocks++;
338 
339  RCP<BlockedCrsMatrix> bA2 = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
341  "Matrix A must be of type BlockedCrsMatrix.");
342 
343  merged2Mat = bA2->Merge();
344 
345  // Add Dirichlet rows to the list of seeds
346  ArrayRCP<const bool> boundaryNodes = Utilities::DetectDirichletRows(*merged2Mat, 0.0);
347  bool haveBoundary = false;
348  for (LO i = 0; i < boundaryNodes.size(); i++)
349  if (boundaryNodes[i]) {
350  // FIXME:
351  // 1. would not this [] overlap with some in the previos blockSeed loop?
352  // 2. do we need to distinguish between pressure and velocity Dirichlet b.c.
353  blockSeeds[i] = numBlocks;
354  haveBoundary = true;
355  }
356  if (haveBoundary)
357  numBlocks++;
358 
359  subList.set("partitioner: type", "user");
360  subList.set("partitioner: map", blockSeeds);
361  subList.set("partitioner: local parts", as<int>(numBlocks));
362 
363  } else {
364  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
365  if (!bA.is_null()) {
366  isBlockedMatrix = true;
367  merged2Mat = bA->Merge();
368  }
369  }
370  }
371 
373  if (isBlockedMatrix == true)
374  tA = Xpetra::toTpetraRowMatrix(merged2Mat);
375  else
376  tA = Xpetra::toTpetraRowMatrix(A_);
377 
378  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
379  SetPrecParameters();
380 
381  prec_->initialize();
382  }
383 
384  prec_->compute();
385 }
386 
387 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
389  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
390 
391  if (this->IsSetup() == true) {
392  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupAggregate(): Setup() has already been called" << std::endl;
393  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupAggregate(): reuse of this type is not available, reverting to full construction" << std::endl;
394  }
395 
396  this->GetOStream(Statistics0) << "Ifpack2Smoother: Using Aggregate Smoothing" << std::endl;
397 
398  RCP<Aggregates> aggregates = Factory::Get<RCP<Aggregates>>(currentLevel, "Aggregates");
399 
400  RCP<const LOMultiVector> vertex2AggId = aggregates->GetVertex2AggId();
401  ArrayRCP<LO> aggregate_ids = rcp_const_cast<LOMultiVector>(vertex2AggId)->getDataNonConst(0);
402  ArrayRCP<LO> dof_ids;
403 
404  // We need to unamalgamate, if the FixedBlockSize > 1
405  LO blocksize = 1;
406  RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
407  if (!matA.is_null())
408  blocksize = matA->GetFixedBlockSize();
409  if (blocksize > 1) {
410  dof_ids.resize(aggregate_ids.size() * blocksize);
411  for (LO i = 0; i < (LO)aggregate_ids.size(); i++) {
412  for (LO j = 0; j < (LO)blocksize; j++)
413  dof_ids[i * blocksize + j] = aggregate_ids[i];
414  }
415  } else {
416  dof_ids = aggregate_ids;
417  }
418 
419  paramList.set("partitioner: map", dof_ids);
420  paramList.set("partitioner: type", "user");
421  paramList.set("partitioner: overlap", 0);
422  paramList.set("partitioner: local parts", (int)aggregates->GetNumAggregates());
423 
424  RCP<const Tpetra::RowMatrix<SC, LO, GO, NO>> tA = Xpetra::toTpetraRowMatrix(A_);
425 
426  type_ = "BLOCKRELAXATION";
427  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
428  SetPrecParameters();
429  prec_->initialize();
430  prec_->compute();
431 }
432 
433 #ifdef HAVE_MUELU_INTREPID2
434 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
436  /*
437 
438  basic notion:
439 
440  Look for user input indicating topo dimension, something like "topological domain type: {node|edge|face}"
441  Call something like what you can find in Poisson example line 1180 to set seeds for a smoother
442 
443  */
444  if (this->IsSetup() == true) {
445  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupTopological(): Setup() has already been called" << std::endl;
446  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupTopological(): reuse of this type is not available, reverting to full construction" << std::endl;
447  }
448 
449  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
450 
451  typedef typename Node::device_type::execution_space ES;
452 
453  typedef Kokkos::DynRankView<LocalOrdinal, typename Node::device_type> FCO; //
454 
456 
457  using namespace std;
458 
459  const Teuchos::RCP<FCO> elemToNode = Factory::Get<Teuchos::RCP<FCO>>(currentLevel, "pcoarsen: element to node map");
460 
461  string basisString = paramList.get<string>("pcoarsen: hi basis");
462  int degree;
463  // NOTE: To make sure Stokhos works we only instantiate these guys with double. There's a lot
464  // of stuff in the guts of Intrepid2 that doesn't play well with Stokhos as of yet. Here, we only
465  // care about the assignment of basis ordinals to topological entities, so this code is actually
466  // independent of the Scalar type--hard-coding double here won't hurt us.
467  auto basis = MueLuIntrepid::BasisFactory<double, ES>(basisString, degree);
468 
469  string topologyTypeString = paramList.get<string>("smoother: neighborhood type");
470  int dimension;
471  if (topologyTypeString == "node")
472  dimension = 0;
473  else if (topologyTypeString == "edge")
474  dimension = 1;
475  else if (topologyTypeString == "face")
476  dimension = 2;
477  else if (topologyTypeString == "cell")
478  dimension = basis->getBaseCellTopology().getDimension();
479  else
480  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unrecognized smoother neighborhood type. Supported types are node, edge, face.");
481  RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_, true);
482  vector<vector<LocalOrdinal>> seeds;
483  MueLuIntrepid::FindGeometricSeedOrdinals(basis, *elemToNode, seeds, *matA->getRowMap(), *matA->getColMap());
484 
485  // Ifpack2 wants the seeds in an array of the same length as the number of local elements,
486  // with local partition #s marked for the ones that are seeds, and invalid for the rest
487  int myNodeCount = matA->getRowMap()->getLocalNumElements();
488  ArrayRCP<LocalOrdinal> nodeSeeds(myNodeCount, lo_invalid);
489  int localPartitionNumber = 0;
490  for (LocalOrdinal seed : seeds[dimension]) {
491  nodeSeeds[seed] = localPartitionNumber++;
492  }
493 
494  paramList.remove("smoother: neighborhood type");
495  paramList.remove("pcoarsen: hi basis");
496 
497  paramList.set("partitioner: map", nodeSeeds);
498  paramList.set("partitioner: type", "user");
499  paramList.set("partitioner: overlap", 1);
500  paramList.set("partitioner: local parts", int(seeds[dimension].size()));
501 
502  RCP<const Tpetra::RowMatrix<SC, LO, GO, NO>> tA = Xpetra::toTpetraRowMatrix(A_);
503 
504  type_ = "BLOCKRELAXATION";
505  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
506  SetPrecParameters();
507  prec_->initialize();
508  prec_->compute();
509 }
510 #endif
511 
512 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
514  if (this->IsSetup() == true) {
515  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupLineSmoothing(): Setup() has already been called" << std::endl;
516  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupLineSmoothing(): reuse of this type is not available, reverting to full construction" << std::endl;
517  }
518 
519  ParameterList& myparamList = const_cast<ParameterList&>(this->GetParameterList());
520 
521  LO CoarseNumZLayers = Factory::Get<LO>(currentLevel, "CoarseNumZLayers");
522  if (CoarseNumZLayers > 0) {
523  Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get<Teuchos::ArrayRCP<LO>>(currentLevel, "LineDetection_VertLineIds");
524 
525  // determine number of local parts
526  LO maxPart = 0;
527  for (size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
528  if (maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
529  }
530  RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_, true);
531  size_t numLocalRows = matA->getLocalNumRows();
532 
533  TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % TVertLineIdSmoo.size() != 0, Exceptions::RuntimeError,
534  "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
535 
536  // actualDofsPerNode is the actual number of matrix rows per mesh element.
537  // It is encoded in either the MueLu Level, or in the Xpetra matrix block size.
538  // This value is needed by Ifpack2 to do decoupled block relaxation.
539  int actualDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
540  LO matrixBlockSize = matA->GetFixedBlockSize();
541  if (matrixBlockSize > 1 && actualDofsPerNode > 1) {
542  TEUCHOS_TEST_FOR_EXCEPTION(actualDofsPerNode != matrixBlockSize, Exceptions::RuntimeError,
543  "MueLu::Ifpack2Smoother::Setup(): A is a block matrix but its block size and DOFs/node from partitioner disagree");
544  } else if (matrixBlockSize > 1) {
545  actualDofsPerNode = matrixBlockSize;
546  }
547  myparamList.set("partitioner: PDE equations", actualDofsPerNode);
548 
549  if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
550  myparamList.set("partitioner: type", "user");
551  myparamList.set("partitioner: map", TVertLineIdSmoo);
552  myparamList.set("partitioner: local parts", maxPart + 1);
553  } else {
554  if (myparamList.isParameter("partitioner: block size") &&
555  myparamList.get<int>("partitioner: block size") != -1) {
556  int block_size = myparamList.get<int>("partitioner: block size");
557  TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % block_size != 0, Exceptions::RuntimeError,
558  "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the specified block size.");
559  numLocalRows /= block_size;
560  }
561 
562  // we assume a constant number of DOFs per node
563  size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
564 
565  // Create a new Teuchos::ArrayRCP<LO> of size numLocalRows and fill it with the corresponding information
567  for (size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
568  for (size_t dof = 0; dof < numDofsPerNode; dof++)
569  partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
570  myparamList.set("partitioner: type", "user");
571  myparamList.set("partitioner: map", partitionerMap);
572  myparamList.set("partitioner: local parts", maxPart + 1);
573  }
574 
575  if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
576  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
577  type_ == "LINESMOOTHING_BANDEDRELAXATION")
578  type_ = "BANDEDRELAXATION";
579  else if (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
580  type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
581  type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
582  type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
583  type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
584  type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION")
585  type_ = "TRIDIAGONALRELAXATION";
586  else
587  type_ = "BLOCKRELAXATION";
588  } else {
589  // line detection failed -> fallback to point-wise relaxation
590  this->GetOStream(Runtime0) << "Line detection failed: fall back to point-wise relaxation" << std::endl;
591  myparamList.remove("partitioner: type", false);
592  myparamList.remove("partitioner: map", false);
593  myparamList.remove("partitioner: local parts", false);
594  type_ = "RELAXATION";
595  }
596 
597  RCP<const Tpetra::RowMatrix<SC, LO, GO, NO>> tA = Xpetra::toTpetraRowMatrix(A_);
598 
599  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
600  SetPrecParameters();
601  prec_->initialize();
602  prec_->compute();
603 }
604 
605 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
607  typedef Tpetra::RowMatrix<SC, LO, GO, NO> tRowMatrix;
608 
609  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
610  if (!bA.is_null())
611  A_ = bA->Merge();
612 
613  RCP<const tRowMatrix> tA = Xpetra::toTpetraRowMatrix(A_);
614 
615  bool reusePreconditioner = false;
616  if (this->IsSetup() == true) {
617  // Reuse the constructed preconditioner
618  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupBlockRelaxation(): Setup() has already been called, assuming reuse" << std::endl;
619 
621  if (!prec.is_null()) {
622  prec->setMatrix(tA);
623  reusePreconditioner = true;
624  } else {
625  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupBlockRelaxation(): reuse of this type is not available (failed cast to CanChangeMatrix), "
626  "reverting to full construction"
627  << std::endl;
628  }
629  }
630 
631  if (!reusePreconditioner) {
632  ParameterList& myparamList = const_cast<ParameterList&>(this->GetParameterList());
633  myparamList.print();
634  if (myparamList.isParameter("partitioner: type") &&
635  myparamList.get<std::string>("partitioner: type") == "line") {
637  Factory::Get<Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO>>>(currentLevel, "Coordinates");
639 
640  RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
641  size_t lclSize = A_->getRangeMap()->getLocalNumElements();
642  if (!matA.is_null())
643  lclSize = matA->getLocalNumRows();
644  size_t numDofsPerNode = lclSize / xCoordinates->getMap()->getLocalNumElements();
645  myparamList.set("partitioner: coordinates", coordinates);
646  myparamList.set("partitioner: PDE equations", (int)numDofsPerNode);
647  }
648 
649  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
650  SetPrecParameters();
651  prec_->initialize();
652  }
653 
654  prec_->compute();
655 }
656 
657 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
659  typedef Tpetra::RowMatrix<SC, LO, GO, NO> tRowMatrix;
660  using STS = Teuchos::ScalarTraits<SC>;
661  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
662  if (!bA.is_null())
663  A_ = bA->Merge();
664 
665  RCP<const tRowMatrix> tA = Xpetra::toTpetraRowMatrix(A_);
666 
667  bool reusePreconditioner = false;
668 
669  if (this->IsSetup() == true) {
670  // Reuse the constructed preconditioner
671  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupChebyshev(): Setup() has already been called, assuming reuse" << std::endl;
672 
674  if (!prec.is_null()) {
675  prec->setMatrix(tA);
676  reusePreconditioner = true;
677  } else {
678  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupChebyshev(): reuse of this type is not available (failed cast to CanChangeMatrix), "
679  "reverting to full construction"
680  << std::endl;
681  }
682  }
683 
684  // Take care of the eigenvalues
685  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
686  SC negone = -STS::one();
687  SC lambdaMax = SetupChebyshevEigenvalues(currentLevel, "A", "", paramList);
688 
689  if (!reusePreconditioner) {
690  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
691  SetPrecParameters();
692  {
693  SubFactoryMonitor m(*this, "Preconditioner init", currentLevel);
694  prec_->initialize();
695  }
696  } else
697  SetPrecParameters();
698 
699  {
700  SubFactoryMonitor m(*this, "Preconditioner compute", currentLevel);
701  prec_->compute();
702  }
703 
704  if (lambdaMax == negone) {
705  typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
706 
708  if (chebyPrec != Teuchos::null) {
709  lambdaMax = chebyPrec->getLambdaMaxForApply();
710  RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
711  if (!matA.is_null())
712  matA->SetMaxEigenvalueEstimate(lambdaMax);
713  this->GetOStream(Statistics1) << "chebyshev: max eigenvalue (calculated by Ifpack2)"
714  << " = " << lambdaMax << std::endl;
715  }
716  TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == negone, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Setup(): no maximum eigenvalue estimate");
717  }
718 }
719 
720 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
722  typedef Tpetra::RowMatrix<SC, LO, GO, NO> tRowMatrix;
723  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
724  if (!bA.is_null())
725  A_ = bA->Merge();
726 
727  RCP<const tRowMatrix> tA = Xpetra::toTpetraRowMatrix(A_);
728 
729  bool reusePreconditioner = false;
730  if (this->IsSetup() == true) {
731  // Reuse the constructed preconditioner
732  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupHiptmair(): Setup() has already been called, assuming reuse" << std::endl;
733 
735  if (!prec.is_null()) {
736  prec->setMatrix(tA);
737  reusePreconditioner = true;
738  } else {
739  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupHiptmair(): reuse of this type is not available (failed cast to CanChangeMatrix), "
740  "reverting to full construction"
741  << std::endl;
742  }
743  }
744 
745  // If we're doing Chebyshev subsmoothing, we'll need to get the eigenvalues
747  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
748  std::string smoother1 = paramList.get("hiptmair: smoother type 1", "CHEBYSHEV");
749  std::string smoother2 = paramList.get("hiptmair: smoother type 2", "CHEBYSHEV");
750  SC lambdaMax11 = negone;
751 
752  if (smoother1 == "CHEBYSHEV") {
753  ParameterList& list1 = paramList.sublist("hiptmair: smoother list 1");
754  lambdaMax11 = SetupChebyshevEigenvalues(currentLevel, "A", "EdgeMatrix ", list1);
755  }
756  if (smoother2 == "CHEBYSHEV") {
757  ParameterList& list2 = paramList.sublist("hiptmair: smoother list 2");
758  SetupChebyshevEigenvalues(currentLevel, "NodeMatrix", "NodeMatrix ", list2);
759  }
760 
761  // FIXME: Should really add some checks to make sure the eigenvalue calcs worked like in
762  // the regular SetupChebyshev
763 
764  // Grab the auxillary matrices and stick them on the list
765  RCP<Operator> NodeMatrix = currentLevel.Get<RCP<Operator>>("NodeMatrix");
766  RCP<Operator> D0 = currentLevel.Get<RCP<Operator>>("D0");
767 
768  RCP<tRowMatrix> tNodeMatrix = Xpetra::toTpetraRowMatrix(NodeMatrix);
769  RCP<tRowMatrix> tD0 = Xpetra::toTpetraRowMatrix(D0);
770 
771  Teuchos::ParameterList newlist;
772  newlist.set("P", tD0);
773  newlist.set("PtAP", tNodeMatrix);
774  if (!reusePreconditioner) {
775  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
776  SetPrecParameters(newlist);
777  prec_->initialize();
778  }
779 
780  prec_->compute();
781 
782  // Post-processing the (1,1) eigenvalue, if we have one
783  if (smoother1 == "CHEBYSHEV" && lambdaMax11 == negone) {
784  using Teuchos::rcp_dynamic_cast;
785  typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
786  auto hiptmairPrec = rcp_dynamic_cast<Ifpack2::Hiptmair<MatrixType>>(prec_);
787  if (hiptmairPrec != Teuchos::null) {
788  auto chebyPrec = rcp_dynamic_cast<Ifpack2::Chebyshev<MatrixType>>(hiptmairPrec->getPrec1());
789  if (chebyPrec != Teuchos::null) {
790  lambdaMax11 = chebyPrec->getLambdaMaxForApply();
791  RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
792  if (!matA.is_null())
793  matA->SetMaxEigenvalueEstimate(lambdaMax11);
794  this->GetOStream(Statistics1) << "chebyshev: max eigenvalue (calculated by Ifpack2)"
795  << " = " << lambdaMax11 << std::endl;
796  }
797  TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax11 == negone, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Setup(): no maximum eigenvalue estimate");
798  }
799  }
800 }
801 
802 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
803 Scalar Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetupChebyshevEigenvalues(Level& currentLevel, const std::string& matrixName, const std::string& label, ParameterList& paramList) const {
804  // Helper: This gets used for smoothers that want to set up Chebyhev
805  typedef Teuchos::ScalarTraits<SC> STS;
806  SC negone = -STS::one();
807  RCP<Operator> currentA = currentLevel.Get<RCP<Operator>>(matrixName);
808  RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(currentA);
809  SC lambdaMax = negone;
810 
811  std::string maxEigString = "chebyshev: max eigenvalue";
812  std::string eigRatioString = "chebyshev: ratio eigenvalue";
813 
814  // Get/calculate the maximum eigenvalue
815  if (paramList.isParameter(maxEigString)) {
816  if (paramList.isType<double>(maxEigString))
817  lambdaMax = paramList.get<double>(maxEigString);
818  else
819  lambdaMax = paramList.get<SC>(maxEigString);
820  this->GetOStream(Statistics1) << label << maxEigString << " (cached with smoother parameter list) = " << lambdaMax << std::endl;
821  if (!matA.is_null())
822  matA->SetMaxEigenvalueEstimate(lambdaMax);
823 
824  } else if (!matA.is_null()) {
825  lambdaMax = matA->GetMaxEigenvalueEstimate();
826  if (lambdaMax != negone) {
827  this->GetOStream(Statistics1) << label << maxEigString << " (cached with matrix) = " << lambdaMax << std::endl;
828  paramList.set(maxEigString, lambdaMax);
829  }
830  }
831 
832  // Calculate the eigenvalue ratio
833  const SC defaultEigRatio = 20;
834 
835  SC ratio = defaultEigRatio;
836  if (paramList.isParameter(eigRatioString)) {
837  if (paramList.isType<double>(eigRatioString))
838  ratio = paramList.get<double>(eigRatioString);
839  else
840  ratio = paramList.get<SC>(eigRatioString);
841  }
842  if (currentLevel.GetLevelID()) {
843  // Update ratio to be
844  // ratio = max(number of fine DOFs / number of coarse DOFs, defaultValue)
845  //
846  // NOTE: We don't need to request previous level matrix as we know for sure it was constructed
847  RCP<const Operator> fineA = currentLevel.GetPreviousLevel()->Get<RCP<Operator>>(matrixName);
848  size_t nRowsFine = fineA->getDomainMap()->getGlobalNumElements();
849  size_t nRowsCoarse = currentA->getDomainMap()->getGlobalNumElements();
850 
851  SC levelRatio = as<SC>(as<double>(nRowsFine) / nRowsCoarse);
852  if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
853  ratio = levelRatio;
854  }
855 
856  this->GetOStream(Statistics1) << label << eigRatioString << " (computed) = " << ratio << std::endl;
857  paramList.set(eigRatioString, ratio);
858 
859  if (paramList.isParameter("chebyshev: use rowsumabs diagonal scaling")) {
860  this->GetOStream(Runtime1) << "chebyshev: using rowsumabs diagonal scaling" << std::endl;
861  bool doScale = false;
862  doScale = paramList.get<bool>("chebyshev: use rowsumabs diagonal scaling");
863  paramList.remove("chebyshev: use rowsumabs diagonal scaling");
864  double chebyReplaceTol = Teuchos::ScalarTraits<Scalar>::eps() * 100;
865  std::string paramName = "chebyshev: rowsumabs diagonal replacement tolerance";
866  if (paramList.isParameter(paramName)) {
867  chebyReplaceTol = paramList.get<double>(paramName);
868  paramList.remove(paramName);
869  }
870  double chebyReplaceVal = Teuchos::ScalarTraits<double>::zero();
871  paramName = "chebyshev: rowsumabs diagonal replacement value";
872  if (paramList.isParameter(paramName)) {
873  chebyReplaceVal = paramList.get<double>(paramName);
874  paramList.remove(paramName);
875  }
876  bool chebyReplaceSingleEntryRowWithZero = false;
877  paramName = "chebyshev: rowsumabs replace single entry row with zero";
878  if (paramList.isParameter(paramName)) {
879  chebyReplaceSingleEntryRowWithZero = paramList.get<bool>(paramName);
880  paramList.remove(paramName);
881  }
882  bool useAverageAbsDiagVal = false;
883  paramName = "chebyshev: rowsumabs use automatic diagonal tolerance";
884  if (paramList.isParameter(paramName)) {
885  useAverageAbsDiagVal = paramList.get<bool>(paramName);
886  paramList.remove(paramName);
887  }
888  if (doScale) {
889  TEUCHOS_ASSERT(!matA.is_null());
890  const bool doReciprocal = true;
891  RCP<Vector> lumpedDiagonal = Utilities::GetLumpedMatrixDiagonal(*matA, doReciprocal, chebyReplaceTol, chebyReplaceVal, chebyReplaceSingleEntryRowWithZero, useAverageAbsDiagVal);
893  paramList.set("chebyshev: operator inv diagonal", tmpVec.getTpetra_Vector());
894  }
895  }
896 
897  return lambdaMax;
898 }
899 
900 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
902  typedef Tpetra::RowMatrix<SC, LO, GO, NO> tRowMatrix;
903  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
904  if (!bA.is_null())
905  A_ = bA->Merge();
906 
907  RCP<const tRowMatrix> tA = Xpetra::toTpetraRowMatrix(A_);
908 
909  bool reusePreconditioner = false;
910  if (this->IsSetup() == true) {
911  // Reuse the constructed preconditioner
912  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
913 
915  if (!prec.is_null()) {
916  prec->setMatrix(tA);
917  reusePreconditioner = true;
918  } else {
919  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupGeneric(): reuse of this type is not available (failed cast to CanChangeMatrix), "
920  "reverting to full construction"
921  << std::endl;
922  }
923  }
924 
925  if (!reusePreconditioner) {
926  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
927  SetPrecParameters();
928  prec_->initialize();
929  }
930 
931  prec_->compute();
932 }
933 
934 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
935 void Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
936  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Apply(): Setup() has not been called");
937 
938  // Forward the InitialGuessIsZero option to Ifpack2
939  // TODO: It might be nice to switch back the internal
940  // "zero starting solution" option of the ifpack2 object prec_ to his
941  // initial value at the end but there is no way right now to get
942  // the current value of the "zero starting solution" in ifpack2.
943  // It's not really an issue, as prec_ can only be used by this method.
944  bool supportInitialGuess = false;
945 
946  if (prec_->supportsZeroStartingSolution()) {
947  prec_->setZeroStartingSolution(InitialGuessIsZero);
948  supportInitialGuess = true;
949  } else if (type_ == "SCHWARZ") {
950  Teuchos::ParameterList paramList;
951  paramList.set("schwarz: zero starting solution", InitialGuessIsZero);
952  // Because additive Schwarz has "delta" semantics, it's sufficient to
953  // toggle only the zero initial guess flag, and not pass in already
954  // set parameters. If we call SetPrecParameters, the subdomain solver
955  // will be destroyed.
956  prec_->setParameters(paramList);
957  supportInitialGuess = true;
958  }
959 
960  // TODO JJH 30Apr2014 Calling SetPrecParameters(paramList) when the smoother
961  // is Ifpack2::AdditiveSchwarz::setParameterList() will destroy the subdomain
962  //(aka inner) solver. This behavior is documented but a departure from what
963  // it previously did, and what other Ifpack2 solvers currently do. So I have
964  // moved SetPrecParameters(paramList) into the if-else block above.
965 
966  // Apply
967  if (InitialGuessIsZero || supportInitialGuess) {
970  prec_->apply(tpB, tpX);
971  } else {
972  using TST = Teuchos::ScalarTraits<Scalar>;
973 
974  RCP<MultiVector> Residual;
975  {
976  std::string prefix = this->ShortClassName() + ": Apply: ";
977  RCP<TimeMonitor> tM = rcp(new TimeMonitor(*this, prefix + "residual calculation", Timings0));
978  Residual = Utilities::Residual(*A_, X, B);
979  }
980 
981  RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
982 
983  Tpetra::MultiVector<SC, LO, GO, NO>& tpX = toTpetra(*Correction);
984  const Tpetra::MultiVector<SC, LO, GO, NO>& tpB = toTpetra(*Residual);
985 
986  prec_->apply(tpB, tpX);
987 
988  X.update(TST::one(), *Correction, TST::one());
989  }
990 }
991 
992 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
994  RCP<Ifpack2Smoother> smoother = rcp(new Ifpack2Smoother(*this));
995  smoother->SetParameterList(this->GetParameterList());
996  return smoother;
997 }
998 
999 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1001  std::ostringstream out;
1003  out << prec_->description();
1004  } else {
1006  out << "{type = " << type_ << "}";
1007  }
1008  return out.str();
1009 }
1010 
1011 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1014 
1015  if (verbLevel & Parameters0)
1016  out0 << "Prec. type: " << type_ << std::endl;
1017 
1018  if (verbLevel & Parameters1) {
1019  out0 << "Parameter list: " << std::endl;
1020  Teuchos::OSTab tab2(out);
1021  out << this->GetParameterList();
1022  out0 << "Overlap: " << overlap_ << std::endl;
1023  }
1024 
1025  if (verbLevel & External)
1026  if (prec_ != Teuchos::null) {
1027  Teuchos::OSTab tab2(out);
1028  out << *prec_ << std::endl;
1029  }
1030 
1031  if (verbLevel & Debug) {
1032  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
1033  << "-" << std::endl
1034  << "RCP<prec_>: " << prec_ << std::endl;
1035  }
1036 }
1037 
1038 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1040  typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
1041  // NOTE: Only works for a subset of Ifpack2's smoothers
1043  if (!pr.is_null()) return pr->getNodeSmootherComplexity();
1044 
1046  if (!pc.is_null()) return pc->getNodeSmootherComplexity();
1047 
1049  if (!pb.is_null()) return pb->getNodeSmootherComplexity();
1050 
1051  RCP<Ifpack2::ILUT<MatrixType>> pi = rcp_dynamic_cast<Ifpack2::ILUT<MatrixType>>(prec_);
1052  if (!pi.is_null()) return pi->getNodeSmootherComplexity();
1053 
1054  RCP<Ifpack2::RILUK<MatrixType>> pk = rcp_dynamic_cast<Ifpack2::RILUK<MatrixType>>(prec_);
1055  if (!pk.is_null()) return pk->getNodeSmootherComplexity();
1056 
1058 }
1059 
1060 } // namespace MueLu
1061 
1062 #endif // HAVE_MUELU_IFPACK2
1063 #endif // MUELU_IFPACK2SMOOTHER_DEF_HPP
void DeclareInput(Level &currentLevel) const
Input.
Important warning messages (one line)
void SetupGeneric(Level &currentLevel)
Exception indicating invalid cast attempted.
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
void FindGeometricSeedOrdinals(Teuchos::RCP< Basis > basis, const LOFieldContainer &elementToNodeMap, std::vector< std::vector< LocalOrdinal > > &seeds, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &columnMap)
KOKKOS_INLINE_FUNCTION LO GetNumAggregates() const
static magnitudeType eps()
Print external lib objects.
GlobalOrdinal GO
T & get(const std::string &name, T def_value)
friend class Ifpack2Smoother
Constructor.
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.
void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
T * get() const
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
size_type size() const
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
std::string type_
ifpack2-specific key phrase that denote smoother type
void SetupBlockRelaxation(Level &currentLevel)
void resize(const size_type n, const T &val=T())
bool isParameter(const std::string &name) const
void SetupSchwarz(Level &currentLevel)
void SetupLineSmoothing(Level &currentLevel)
Print statistics that do not involve significant additional computation.
bool remove(std::string const &name, bool throwIfNotExists=true)
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
RCP< SmootherPrototype > Copy() const
void declareConstructionOutcome(bool fail, std::string msg)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
bool isSublist(const std::string &name) const
void Setup(Level &currentLevel)
Set up the smoother.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Magnitude >::zero(), bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
void SetupHiptmair(Level &currentLevel)
void SetupTopological(Level &currentLevel)
std::string description() const
Return a simple one-line description of this object.
bool IsSetup() const
Get the state of a smoother prototype.
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)
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
ParameterList & setParameters(const ParameterList &source)
Print all timing information.
Print class parameters.
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
RCP< Level > & GetPreviousLevel()
Previous level.
Definition: MueLu_Level.cpp:60
MatrixType::scalar_type getLambdaMaxForApply() const
Scalar SC
void SetupAggregate(Level &currentLevel)
bool isType(const std::string &name) const
RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_Vector() const
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
Node NO
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos::RCP< BlockCrsMatrix< Scalar, LO, GO, Node > > convertToBlockCrsMatrix(const Tpetra::CrsMatrix< Scalar, LO, GO, Node > &pointMatrix, const LO &blockSize, bool use_local_ID=true)
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:51
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
void SetupChebyshev(Level &currentLevel)
#define TEUCHOS_ASSERT(assertion_test)
Description of what is happening (more verbose)
Class that encapsulates Ifpack2 smoothers.
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
virtual std::string description() const
Return a simple one-line description of this object.
std::string toString(const T &t)
bool is_null() const
Scalar SetupChebyshevEigenvalues(Level &currentLevel, const std::string &matrixName, const std::string &label, Teuchos::ParameterList &paramList) const