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_Relaxation.hpp>
59 #include <Ifpack2_ILUT.hpp>
60 #include <Ifpack2_BlockRelaxation.hpp>
61 #include <Ifpack2_Factory.hpp>
62 #include <Ifpack2_Parameters.hpp>
63 
64 #include <Xpetra_BlockedCrsMatrix.hpp>
65 #include <Xpetra_CrsMatrix.hpp>
66 #include <Xpetra_CrsMatrixWrap.hpp>
67 #include <Xpetra_Matrix.hpp>
68 #include <Xpetra_MultiVectorFactory.hpp>
69 #include <Xpetra_TpetraMultiVector.hpp>
70 
72 #include "MueLu_Level.hpp"
74 #include "MueLu_Utilities.hpp"
75 #include "MueLu_Monitor.hpp"
76 
77 #ifdef HAVE_MUELU_INTREPID2
80 #include "Intrepid2_Basis.hpp"
81 #include "Kokkos_DynRankView.hpp"
82 #endif
83 
84 // #define IFPACK2_HAS_PROPER_REUSE
85 
86 namespace MueLu {
87 
88  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
89  Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Ifpack2Smoother(const std::string& type, const Teuchos::ParameterList& paramList, const LO& overlap)
90  : type_(type), overlap_(overlap)
91  {
92  SetParameterList(paramList);
93  }
94 
95  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
97  Factory::SetParameterList(paramList);
98 
100  // It might be invalid to change parameters after the setup, but it depends entirely on Ifpack implementation.
101  // TODO: I don't know if Ifpack returns an error code or exception or ignore parameters modification in this case...
102  SetPrecParameters();
103  }
104  }
105 
106  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
108  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
109  paramList.setParameters(list);
110 
111  RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
112 
113  prec_->setParameters(*precList);
114 
115  paramList.setParameters(*precList); // what about that??
116  }
117 
118  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
120  this->Input(currentLevel, "A");
121 
122  if (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
123  type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
124  type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
125  type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
126  type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
127  type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION" ||
128  type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
129  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
130  type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
131  type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
132  type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
133  type_ == "LINESMOOTHING_BLOCKRELAXATION") {
134  this->Input(currentLevel, "CoarseNumZLayers"); // necessary for fallback criterion
135  this->Input(currentLevel, "LineDetection_VertLineIds"); // necessary to feed block smoother
136  }
137  else if (type_ == "BLOCK RELAXATION" ||
138  type_ == "BLOCK_RELAXATION" ||
139  type_ == "BLOCKRELAXATION")
140  {
141  //We need to check for the "partitioner type" = "line"
142  ParameterList precList = this->GetParameterList();
143  if(precList.isParameter("partitioner: type") &&
144  precList.get<std::string>("partitioner: type") == "line") {
145  this->Input(currentLevel, "Coordinates");
146  }
147  }
148  else if (type_ == "TOPOLOGICAL")
149  {
150  // for the topological smoother, we require an element to node map:
151  this->Input(currentLevel, "pcoarsen: element to node map");
152  }
153  }
154 
155  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
157  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
158 
159  A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");
160 
161  if (type_ == "SCHWARZ")
162  SetupSchwarz(currentLevel);
163 
164  else if (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
165  type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
166  type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
167  type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
168  type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
169  type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION" ||
170  type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
171  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
172  type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
173  type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
174  type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
175  type_ == "LINESMOOTHING_BLOCKRELAXATION")
176  SetupLineSmoothing(currentLevel);
177 
178  else if (type_ == "BLOCK_RELAXATION" ||
179  type_ == "BLOCK RELAXATION" ||
180  type_ == "BLOCKRELAXATION")
181  SetupBlockRelaxation(currentLevel);
182 
183  else if (type_ == "CHEBYSHEV")
184  SetupChebyshev(currentLevel);
185 
186  else if (type_ == "TOPOLOGICAL")
187  {
188 #ifdef HAVE_MUELU_INTREPID2
189  SetupTopological(currentLevel);
190 #else
191  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "'TOPOLOGICAL' smoother choice requires Intrepid2");
192 #endif
193  }
194  else
195  {
196  SetupGeneric(currentLevel);
197  }
198 
200 
201  this->GetOStream(Statistics1) << description() << std::endl;
202  }
203 
204  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
206  typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
207 
208  bool reusePreconditioner = false;
209  if (this->IsSetup() == true) {
210  // Reuse the constructed preconditioner
211  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupSchwarz(): Setup() has already been called, assuming reuse" << std::endl;
212 
213  bool isTRowMatrix = true;
214  RCP<const tRowMatrix> tA;
215  try {
217  } catch (Exceptions::BadCast) {
218  isTRowMatrix = false;
219  }
220 
221  RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
222  if (!prec.is_null() && isTRowMatrix) {
223 #ifdef IFPACK2_HAS_PROPER_REUSE
224  prec->resetMatrix(tA);
225  reusePreconditioner = true;
226 #else
227  this->GetOStream(Errors) << "Ifpack2 does not have proper reuse yet." << std::endl;
228 #endif
229 
230  } else {
231  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available "
232  "(either failed cast to CanChangeMatrix, or to Tpetra Row Matrix), reverting to full construction" << std::endl;
233  }
234  }
235 
236  if (!reusePreconditioner) {
237  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
238 
239  bool isBlockedMatrix = false;
240  RCP<Matrix> merged2Mat;
241 
242  std::string sublistName = "subdomain solver parameters";
243  if (paramList.isSublist(sublistName)) {
244  // If we are doing "user" partitioning, we assume that what the user
245  // really wants to do is make tiny little subdomains with one row
246  // assigned to each subdomain. The rows used for these little
247  // subdomains correspond to those in the 2nd block row. Then,
248  // if we overlap these mini-subdomains, we will do something that
249  // looks like Vanka (grabbing all velocities associated with each
250  // each pressure unknown). In addition, we put all Dirichlet points
251  // as a little mini-domain.
252  ParameterList& subList = paramList.sublist(sublistName);
253 
254  std::string partName = "partitioner: type";
255  if (subList.isParameter(partName) && subList.get<std::string>(partName) == "user") {
256  isBlockedMatrix = true;
257 
258  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
259  TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
260  "Matrix A must be of type BlockedCrsMatrix.");
261 
262  size_t numVels = bA->getMatrix(0,0)->getNodeNumRows();
263  size_t numPres = bA->getMatrix(1,0)->getNodeNumRows();
264  size_t numRows = A_->getNodeNumRows();
265 
266  ArrayRCP<LocalOrdinal> blockSeeds(numRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
267 
268  size_t numBlocks = 0;
269  for (size_t rowOfB = numVels; rowOfB < numVels+numPres; ++rowOfB)
270  blockSeeds[rowOfB] = numBlocks++;
271 
272  RCP<BlockedCrsMatrix> bA2 = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
273  TEUCHOS_TEST_FOR_EXCEPTION(bA2.is_null(), Exceptions::BadCast,
274  "Matrix A must be of type BlockedCrsMatrix.");
275 
276  merged2Mat = bA2->Merge();
277 
278  // Add Dirichlet rows to the list of seeds
279  ArrayRCP<const bool> boundaryNodes = Utilities::DetectDirichletRows(*merged2Mat, 0.0);
280  bool haveBoundary = false;
281  for (LO i = 0; i < boundaryNodes.size(); i++)
282  if (boundaryNodes[i]) {
283  // FIXME:
284  // 1. would not this [] overlap with some in the previos blockSeed loop?
285  // 2. do we need to distinguish between pressure and velocity Dirichlet b.c.
286  blockSeeds[i] = numBlocks;
287  haveBoundary = true;
288  }
289  if (haveBoundary)
290  numBlocks++;
291 
292  subList.set("partitioner: map", blockSeeds);
293  subList.set("partitioner: local parts", as<int>(numBlocks));
294 
295  } else {
296  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
297  if (!bA.is_null()) {
298  isBlockedMatrix = true;
299  merged2Mat = bA->Merge();
300  }
301  }
302  }
303 
304  RCP<const tRowMatrix> tA;
305  if (isBlockedMatrix == true) tA = Utilities::Op2NonConstTpetraRow(merged2Mat);
306  else tA = Utilities::Op2NonConstTpetraRow(A_);
307 
308  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
309  SetPrecParameters();
310 
311  prec_->initialize();
312  }
313 
314  prec_->compute();
315  }
316 
317 #ifdef HAVE_MUELU_INTREPID2
318  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
320  /*
321 
322  basic notion:
323 
324  Look for user input indicating topo dimension, something like "topological domain type: {node|edge|face}"
325  Call something like what you can find in Poisson example line 1180 to set seeds for a smoother
326 
327  */
328  if (this->IsSetup() == true) {
329  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupTopological(): Setup() has already been called" << std::endl;
330  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupTopological(): reuse of this type is not available, reverting to full construction" << std::endl;
331  }
332 
333  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
334 
335  typedef typename Node::device_type::execution_space ES;
336 
337  typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCO; //
338 
339  LocalOrdinal lo_invalid = Teuchos::OrdinalTraits<LO>::invalid();
340 
341  using namespace std;
342 
343  const Teuchos::RCP<FCO> elemToNode = Factory::Get<Teuchos::RCP<FCO> >(currentLevel,"pcoarsen: element to node map");
344 
345  string basisString = paramList.get<string>("pcoarsen: hi basis");
346  int degree;
347  // NOTE: To make sure Stokhos works we only instantiate these guys with double. There's a lot
348  // of stuff in the guts of Intrepid2 that doesn't play well with Stokhos as of yet. Here, we only
349  // care about the assignment of basis ordinals to topological entities, so this code is actually
350  // independent of the Scalar type--hard-coding double here won't hurt us.
351  auto basis = MueLuIntrepid::BasisFactory<double,ES>(basisString, degree);
352 
353  string topologyTypeString = paramList.get<string>("smoother: neighborhood type");
354  int dimension;
355  if (topologyTypeString == "node")
356  dimension = 0;
357  else if (topologyTypeString == "edge")
358  dimension = 1;
359  else if (topologyTypeString == "face")
360  dimension = 2;
361  else if (topologyTypeString == "cell")
362  dimension = basis->getBaseCellTopology().getDimension();
363  else
364  TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unrecognized smoother neighborhood type. Supported types are node, edge, face.");
365  vector<vector<LocalOrdinal>> seeds;
366  MueLuIntrepid::FindGeometricSeedOrdinals(basis, *elemToNode, seeds, *A_->getRowMap(), *A_->getColMap());
367 
368  // Ifpack2 wants the seeds in an array of the same length as the number of local elements,
369  // with local partition #s marked for the ones that are seeds, and invalid for the rest
370  int myNodeCount = A_->getRowMap()->getNodeNumElements();
371  ArrayRCP<LocalOrdinal> nodeSeeds(myNodeCount,lo_invalid);
372  int localPartitionNumber = 0;
373  for (LocalOrdinal seed : seeds[dimension])
374  {
375  nodeSeeds[seed] = localPartitionNumber++;
376  }
377 
378  paramList.remove("smoother: neighborhood type");
379  paramList.remove("pcoarsen: hi basis");
380 
381  paramList.set("partitioner: map", nodeSeeds);
382  paramList.set("partitioner: type", "user");
383  paramList.set("partitioner: overlap", 1);
384  paramList.set("partitioner: local parts", int(seeds[dimension].size()));
385 
386  RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tA = Utilities::Op2NonConstTpetraRow(A_);
387 
388  type_ = "BLOCKRELAXATION";
389  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
390  SetPrecParameters();
391  prec_->initialize();
392  prec_->compute();
393  }
394 #endif
395 
396  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
398  if (this->IsSetup() == true) {
399  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupLineSmoothing(): Setup() has already been called" << std::endl;
400  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupLineSmoothing(): reuse of this type is not available, reverting to full construction" << std::endl;
401  }
402 
403  ParameterList& myparamList = const_cast<ParameterList&>(this->GetParameterList());
404 
405  LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,"CoarseNumZLayers");
406  if (CoarseNumZLayers > 0) {
407  Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get< Teuchos::ArrayRCP<LO> >(currentLevel, "LineDetection_VertLineIds");
408 
409  // determine number of local parts
410  LO maxPart = 0;
411  for(size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
412  if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
413  }
414 
415  size_t numLocalRows = A_->getNodeNumRows();
416  TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % TVertLineIdSmoo.size() != 0, Exceptions::RuntimeError,
417  "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
418 
419  if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
420  myparamList.set("partitioner: type","user");
421  myparamList.set("partitioner: map",TVertLineIdSmoo);
422  myparamList.set("partitioner: local parts",maxPart+1);
423  } else {
424  // we assume a constant number of DOFs per node
425  size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
426 
427  // Create a new Teuchos::ArrayRCP<LO> of size numLocalRows and fill it with the corresponding information
428  Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
429  for (size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
430  for (size_t dof = 0; dof < numDofsPerNode; dof++)
431  partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
432  myparamList.set("partitioner: type","user");
433  myparamList.set("partitioner: map",partitionerMap);
434  myparamList.set("partitioner: local parts",maxPart + 1);
435  }
436 
437  if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
438  type_ == "LINESMOOTHING_BANDED RELAXATION" ||
439  type_ == "LINESMOOTHING_BANDEDRELAXATION")
440  type_ = "BANDEDRELAXATION";
441  else if (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
442  type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
443  type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
444  type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
445  type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
446  type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION")
447  type_ = "TRIDIAGONALRELAXATION";
448  else
449  type_ = "BLOCKRELAXATION";
450  } else {
451  // line detection failed -> fallback to point-wise relaxation
452  this->GetOStream(Runtime0) << "Line detection failed: fall back to point-wise relaxation" << std::endl;
453  myparamList.remove("partitioner: type",false);
454  myparamList.remove("partitioner: map", false);
455  myparamList.remove("partitioner: local parts",false);
456  type_ = "RELAXATION";
457  }
458 
459  RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tA = Utilities::Op2NonConstTpetraRow(A_);
460 
461  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
462  SetPrecParameters();
463  prec_->initialize();
464  prec_->compute();
465  }
466 
467  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
469  typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
470 
471  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
472  if (!bA.is_null())
473  A_ = bA->Merge();
474 
475  RCP<const tRowMatrix> tA = Utilities::Op2NonConstTpetraRow(A_);
476 
477  bool reusePreconditioner = false;
478  if (this->IsSetup() == true) {
479  // Reuse the constructed preconditioner
480  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
481 
482  RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
483  if (!prec.is_null()) {
484 #ifdef IFPACK2_HAS_PROPER_REUSE
485  prec->resetMatrix(tA);
486  reusePreconditioner = true;
487 #else
488  this->GetOStream(Errors) << "Ifpack2 does not have proper reuse yet." << std::endl;
489 #endif
490 
491  } else {
492  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available (failed cast to CanChangeMatrix), "
493  "reverting to full construction" << std::endl;
494  }
495  }
496 
497  if (!reusePreconditioner) {
498  ParameterList& myparamList = const_cast<ParameterList&>(this->GetParameterList());
499  myparamList.print();
500  if(myparamList.isParameter("partitioner: type") &&
501  myparamList.get<std::string>("partitioner: type") == "line") {
502  Teuchos::RCP<Xpetra::MultiVector<double,LO,GO,NO> > xCoordinates =
503  Factory::Get<Teuchos::RCP<Xpetra::MultiVector<double,LO,GO,NO> > >(currentLevel, "Coordinates");
504  Teuchos::RCP<Tpetra::MultiVector<double,LO,GO,NO> > coordinates = Teuchos::rcpFromRef(Xpetra::toTpetra<double,LO,GO,NO>(*xCoordinates));
505  myparamList.set("partitioner: coordinates", coordinates);
506  }
507 
508  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
509  SetPrecParameters();
510  prec_->initialize();
511  }
512 
513  prec_->compute();
514  }
515 
516  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
518  if (this->IsSetup() == true) {
519  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupChebyshev(): SetupChebyshev() has already been called" << std::endl;
520  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupChebyshev(): reuse of this type is not available, reverting to full construction" << std::endl;
521  }
522 
523  typedef Teuchos::ScalarTraits<SC> STS;
524  SC negone = -STS::one();
525 
526  SC lambdaMax = negone;
527  {
528  std::string maxEigString = "chebyshev: max eigenvalue";
529  std::string eigRatioString = "chebyshev: ratio eigenvalue";
530 
531  ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
532 
533  // Get/calculate the maximum eigenvalue
534  if (paramList.isParameter(maxEigString)) {
535  if (paramList.isType<double>(maxEigString))
536  lambdaMax = paramList.get<double>(maxEigString);
537  else
538  lambdaMax = paramList.get<SC>(maxEigString);
539  this->GetOStream(Statistics1) << maxEigString << " (cached with smoother parameter list) = " << lambdaMax << std::endl;
540 
541  } else {
542  lambdaMax = A_->GetMaxEigenvalueEstimate();
543  if (lambdaMax != negone) {
544  this->GetOStream(Statistics1) << maxEigString << " (cached with matrix) = " << lambdaMax << std::endl;
545  paramList.set(maxEigString, lambdaMax);
546  }
547  }
548 
549  // Calculate the eigenvalue ratio
550  const SC defaultEigRatio = 20;
551 
552  SC ratio = defaultEigRatio;
553  if (paramList.isParameter(eigRatioString)) {
554  if (paramList.isType<double>(eigRatioString))
555  ratio = paramList.get<double>(eigRatioString);
556  else
557  ratio = paramList.get<SC>(eigRatioString);
558  }
559  if (currentLevel.GetLevelID()) {
560  // Update ratio to be
561  // ratio = max(number of fine DOFs / number of coarse DOFs, defaultValue)
562  //
563  // NOTE: We don't need to request previous level matrix as we know for sure it was constructed
564  RCP<const Matrix> fineA = currentLevel.GetPreviousLevel()->Get<RCP<Matrix> >("A");
565  size_t nRowsFine = fineA->getGlobalNumRows();
566  size_t nRowsCoarse = A_->getGlobalNumRows();
567 
568  SC levelRatio = as<SC>(as<float>(nRowsFine)/nRowsCoarse);
569  if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
570  ratio = levelRatio;
571  }
572 
573  this->GetOStream(Statistics1) << eigRatioString << " (computed) = " << ratio << std::endl;
574  paramList.set(eigRatioString, ratio);
575  }
576 
577  RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tA = Utilities::Op2NonConstTpetraRow(A_);
578 
579  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
580  SetPrecParameters();
581  {
582  SubFactoryMonitor(*this, "Preconditioner init", currentLevel);
583  prec_->initialize();
584  }
585  {
586  SubFactoryMonitor(*this, "Preconditioner compute", currentLevel);
587  prec_->compute();
588  }
589 
590  if (lambdaMax == negone) {
591  typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
592 
593  Teuchos::RCP<Ifpack2::Chebyshev<MatrixType> > chebyPrec = rcp_dynamic_cast<Ifpack2::Chebyshev<MatrixType> >(prec_);
594  if (chebyPrec != Teuchos::null) {
595  lambdaMax = chebyPrec->getLambdaMaxForApply();
596  A_->SetMaxEigenvalueEstimate(lambdaMax);
597  this->GetOStream(Statistics1) << "chebyshev: max eigenvalue (calculated by Ifpack2)" << " = " << lambdaMax << std::endl;
598  }
599  TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == negone, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Setup(): no maximum eigenvalue estimate");
600  }
601  }
602 
603  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
605  typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
606 
607  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
608  if (!bA.is_null())
609  A_ = bA->Merge();
610 
611  RCP<const tRowMatrix> tA = Utilities::Op2NonConstTpetraRow(A_);
612 
613  bool reusePreconditioner = false;
614  if (this->IsSetup() == true) {
615  // Reuse the constructed preconditioner
616  this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
617 
618  RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
619  if (!prec.is_null()) {
620 #ifdef IFPACK2_HAS_PROPER_REUSE
621  prec->resetMatrix(tA);
622  reusePreconditioner = true;
623 #else
624  this->GetOStream(Errors) << "Ifpack2 does not have proper reuse yet." << std::endl;
625 #endif
626 
627  } else {
628  this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available (failed cast to CanChangeMatrix), "
629  "reverting to full construction" << std::endl;
630  }
631  }
632 
633  if (!reusePreconditioner) {
634  prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
635  SetPrecParameters();
636  prec_->initialize();
637  }
638 
639  prec_->compute();
640  }
641 
642  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
643  void Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
644  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Apply(): Setup() has not been called");
645 
646  // Forward the InitialGuessIsZero option to Ifpack2
647  // TODO: It might be nice to switch back the internal
648  // "zero starting solution" option of the ifpack2 object prec_ to his
649  // initial value at the end but there is no way right now to get
650  // the current value of the "zero starting solution" in ifpack2.
651  // It's not really an issue, as prec_ can only be used by this method.
652  // TODO: When https://software.sandia.gov/bugzilla/show_bug.cgi?id=5283#c2 is done
653  // we should remove the if/else/elseif and just test if this
654  // option is supported by current ifpack2 preconditioner
655  Teuchos::ParameterList paramList;
656  bool supportInitialGuess = false;
657  if (type_ == "CHEBYSHEV") {
658  paramList.set("chebyshev: zero starting solution", InitialGuessIsZero);
659  SetPrecParameters(paramList);
660  supportInitialGuess = true;
661 
662  } else if (type_ == "RELAXATION") {
663  paramList.set("relaxation: zero starting solution", InitialGuessIsZero);
664  SetPrecParameters(paramList);
665  supportInitialGuess = true;
666 
667  } else if (type_ == "KRYLOV") {
668  paramList.set("krylov: zero starting solution", InitialGuessIsZero);
669  SetPrecParameters(paramList);
670  supportInitialGuess = true;
671 
672  } else if (type_ == "SCHWARZ") {
673  paramList.set("schwarz: zero starting solution", InitialGuessIsZero);
674  //Because additive Schwarz has "delta" semantics, it's sufficient to
675  //toggle only the zero initial guess flag, and not pass in already
676  //set parameters. If we call SetPrecParameters, the subdomain solver
677  //will be destroyed.
678  prec_->setParameters(paramList);
679  supportInitialGuess = true;
680  }
681 
682  //TODO JJH 30Apr2014 Calling SetPrecParameters(paramList) when the smoother
683  //is Ifpack2::AdditiveSchwarz::setParameterList() will destroy the subdomain
684  //(aka inner) solver. This behavior is documented but a departure from what
685  //it previously did, and what other Ifpack2 solvers currently do. So I have
686  //moved SetPrecParameters(paramList) into the if-else block above.
687 
688  // Apply
689  if (InitialGuessIsZero || supportInitialGuess) {
690  Tpetra::MultiVector<SC,LO,GO,NO>& tpX = Utilities::MV2NonConstTpetraMV(X);
691  const Tpetra::MultiVector<SC,LO,GO,NO>& tpB = Utilities::MV2TpetraMV(B);
692  prec_->apply(tpB, tpX);
693  } else {
694  typedef Teuchos::ScalarTraits<Scalar> TST;
695  RCP<MultiVector> Residual = Utilities::Residual(*A_, X, B);
696  RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
697 
698  Tpetra::MultiVector<SC,LO,GO,NO>& tpX = Utilities::MV2NonConstTpetraMV(*Correction);
699  const Tpetra::MultiVector<SC,LO,GO,NO>& tpB = Utilities::MV2TpetraMV(*Residual);
700 
701  prec_->apply(tpB, tpX);
702 
703  X.update(TST::one(), *Correction, TST::one());
704  }
705  }
706 
707  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
708  RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Copy() const {
709  RCP<Ifpack2Smoother> smoother = rcp(new Ifpack2Smoother(*this) );
710  smoother->SetParameterList(this->GetParameterList());
711  return smoother;
712  }
713 
714  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
716  std::ostringstream out;
718  out << prec_->description();
719  } else {
721  out << "{type = " << type_ << "}";
722  }
723  return out.str();
724  }
725 
726  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
727  void Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
729 
730  if (verbLevel & Parameters0)
731  out0 << "Prec. type: " << type_ << std::endl;
732 
733  if (verbLevel & Parameters1) {
734  out0 << "Parameter list: " << std::endl;
735  Teuchos::OSTab tab2(out);
736  out << this->GetParameterList();
737  out0 << "Overlap: " << overlap_ << std::endl;
738  }
739 
740  if (verbLevel & External)
741  if (prec_ != Teuchos::null) {
742  Teuchos::OSTab tab2(out);
743  out << *prec_ << std::endl;
744  }
745 
746  if (verbLevel & Debug) {
747  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
748  << "-" << std::endl
749  << "RCP<prec_>: " << prec_ << std::endl;
750  }
751  }
752 
753  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
755  typedef Tpetra::RowMatrix<SC,LO,GO,NO> MatrixType;
756  // NOTE: Only works for a subset of Ifpack2's smoothers
757  RCP<Ifpack2::Relaxation<MatrixType> > pr = rcp_dynamic_cast<Ifpack2::Relaxation<MatrixType> >(prec_);
758  if(!pr.is_null()) return pr->getNodeSmootherComplexity();
759 
760  RCP<Ifpack2::Chebyshev<MatrixType> > pc = rcp_dynamic_cast<Ifpack2::Chebyshev<MatrixType> >(prec_);
761  if(!pc.is_null()) return pc->getNodeSmootherComplexity();
762 
763  RCP<Ifpack2::BlockRelaxation<MatrixType> > pb = rcp_dynamic_cast<Ifpack2::BlockRelaxation<MatrixType> >(prec_);
764  if(!pb.is_null()) return pb->getNodeSmootherComplexity();
765 
766  RCP<Ifpack2::ILUT<MatrixType> > pi = rcp_dynamic_cast<Ifpack2::ILUT<MatrixType> >(prec_);
767  if(!pi.is_null()) return pi->getNodeSmootherComplexity();
768 
769  RCP<Ifpack2::RILUK<MatrixType> > pk = rcp_dynamic_cast<Ifpack2::RILUK<MatrixType> >(prec_);
770  if(!pk.is_null()) return pk->getNodeSmootherComplexity();
771 
772 
773  return Teuchos::OrdinalTraits<size_t>::invalid();
774  }
775 
776 
777 } // namespace MueLu
778 
779 #endif // HAVE_MUELU_TPETRA && HAVE_MUELU_IFPACK2
780 #endif // MUELU_IFPACK2SMOOTHER_DEF_HPP
Important warning messages (one line)
void SetupGeneric(Level &currentLevel)
RCP< SmootherPrototype > Copy() const
Exception indicating invalid cast attempted.
RCP< Level > & GetPreviousLevel()
Previous level.
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
size_t getNodeSmootherComplexity() const
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)
Namespace for MueLu classes and methods.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
void SetupBlockRelaxation(Level &currentLevel)
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void SetupSchwarz(Level &currentLevel)
void SetupLineSmoothing(Level &currentLevel)
MatrixType::scalar_type getLambdaMaxForApply() const
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
virtual void SetParameterList(const ParameterList &paramList)
Set parameters from a parameter list and return with default values.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Setup(Level &currentLevel)
Set up the smoother.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
size_t getNodeSmootherComplexity() const
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.
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.