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