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