MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_ShiftedLaplacian_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_SHIFTEDLAPLACIAN_DEF_HPP
11 #define MUELU_SHIFTEDLAPLACIAN_DEF_HPP
12 
13 #include "MueLu_ShiftedLaplacian_decl.hpp"
14 
15 #if defined(HAVE_MUELU_IFPACK2) and defined(HAVE_MUELU_TPETRA)
16 
17 #include <MueLu_AmalgamationFactory.hpp>
18 #include <MueLu_CoalesceDropFactory.hpp>
19 #include <MueLu_CoarseMapFactory.hpp>
20 #include <MueLu_CoupledRBMFactory.hpp>
21 #include <MueLu_DirectSolver.hpp>
22 #include <MueLu_GenericRFactory.hpp>
23 #include <MueLu_Hierarchy.hpp>
24 #include <MueLu_Ifpack2Smoother.hpp>
25 #include <MueLu_PFactory.hpp>
26 #include <MueLu_PgPFactory.hpp>
27 #include <MueLu_RAPFactory.hpp>
28 #include <MueLu_RAPShiftFactory.hpp>
29 #include <MueLu_SaPFactory.hpp>
30 #include <MueLu_ShiftedLaplacian.hpp>
31 #include <MueLu_ShiftedLaplacianOperator.hpp>
32 #include <MueLu_SmootherFactory.hpp>
33 #include <MueLu_SmootherPrototype.hpp>
34 #include <MueLu_TentativePFactory.hpp>
35 #include <MueLu_TransPFactory.hpp>
36 #include <MueLu_UncoupledAggregationFactory.hpp>
37 #include <MueLu_Utilities.hpp>
38 
39 namespace MueLu {
40 
41 // Destructor
42 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
43 ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~ShiftedLaplacian() {}
44 
45 // Input
46 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
47 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameters(Teuchos::RCP<Teuchos::ParameterList> paramList) {
48  // Parameters
49  coarseGridSize_ = paramList->get("MueLu: coarse size", 1000);
50  numLevels_ = paramList->get("MueLu: levels", 3);
51  int stype = paramList->get("MueLu: smoother", 8);
52  if (stype == 1) {
53  Smoother_ = "jacobi";
54  } else if (stype == 2) {
55  Smoother_ = "gauss-seidel";
56  } else if (stype == 3) {
57  Smoother_ = "symmetric gauss-seidel";
58  } else if (stype == 4) {
59  Smoother_ = "chebyshev";
60  } else if (stype == 5) {
61  Smoother_ = "krylov";
62  } else if (stype == 6) {
63  Smoother_ = "ilut";
64  } else if (stype == 7) {
65  Smoother_ = "riluk";
66  } else if (stype == 8) {
67  Smoother_ = "schwarz";
68  } else if (stype == 9) {
69  Smoother_ = "superilu";
70  } else if (stype == 10) {
71  Smoother_ = "superlu";
72  } else {
73  Smoother_ = "schwarz";
74  }
75  smoother_sweeps_ = paramList->get("MueLu: sweeps", 5);
76  smoother_damping_ = paramList->get("MueLu: relax val", 1.0);
77  ncycles_ = paramList->get("MueLu: cycles", 1);
78  iters_ = paramList->get("MueLu: iterations", 500);
79  solverType_ = paramList->get("MueLu: solver type", 1);
80  restart_size_ = paramList->get("MueLu: restart size", 100);
81  recycle_size_ = paramList->get("MueLu: recycle size", 25);
82  isSymmetric_ = paramList->get("MueLu: symmetric", true);
83  ilu_leveloffill_ = paramList->get("MueLu: level-of-fill", 5);
84  ilu_abs_thresh_ = paramList->get("MueLu: abs thresh", 0.0);
85  ilu_rel_thresh_ = paramList->get("MueLu: rel thresh", 1.0);
86  ilu_diagpivotthresh_ = paramList->get("MueLu: piv thresh", 0.1);
87  ilu_drop_tol_ = paramList->get("MueLu: drop tol", 0.01);
88  ilu_fill_tol_ = paramList->get("MueLu: fill tol", 0.01);
89  schwarz_overlap_ = paramList->get("MueLu: overlap", 0);
90  schwarz_usereorder_ = paramList->get("MueLu: use reorder", true);
91  int combinemode = paramList->get("MueLu: combine mode", 1);
92  if (combinemode == 0) {
93  schwarz_combinemode_ = Tpetra::ZERO;
94  } else {
95  schwarz_combinemode_ = Tpetra::ADD;
96  }
97  tol_ = paramList->get("MueLu: tolerance", 0.001);
98 }
99 
100 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
101 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setProblemMatrix(RCP<Matrix>& A) {
102  A_ = A;
103  if (A_ != Teuchos::null)
104  TpetraA_ = Utilities::Op2NonConstTpetraCrs(A_);
105 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
106  if (LinearProblem_ != Teuchos::null)
107  LinearProblem_->setOperator(TpetraA_);
108 #else
109  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "ShiftedLaplacian only available with Tpetra and GO=int enabled.");
110 #endif
111 }
112 
113 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
114 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setProblemMatrix(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraA) {
115  TpetraA_ = TpetraA;
116 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
117  if (LinearProblem_ != Teuchos::null)
118  LinearProblem_->setOperator(TpetraA_);
119 #endif
120 }
121 
122 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
123 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setPreconditioningMatrix(RCP<Matrix>& P) {
124  P_ = P;
125  GridTransfersExist_ = false;
126 }
127 
128 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
129 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setPreconditioningMatrix(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraP) {
130  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Atmp = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraP));
131  P_ = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp));
132  GridTransfersExist_ = false;
133 }
134 
135 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
136 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setstiff(RCP<Matrix>& K) {
137  K_ = K;
138 }
139 
140 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setstiff(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraK) {
142  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Atmp = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraK));
143  K_ = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp));
144 }
145 
146 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
147 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setmass(RCP<Matrix>& M) {
148  M_ = M;
149 }
150 
151 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
152 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setmass(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraM) {
153  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Atmp = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraM));
154  M_ = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp));
155 }
156 
157 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
158 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setcoords(RCP<MultiVector>& Coords) {
159  Coords_ = Coords;
160 }
161 
162 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
163 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setNullSpace(RCP<MultiVector> NullSpace) {
164  NullSpace_ = NullSpace;
165 }
166 
167 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
168 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setLevelShifts(std::vector<Scalar> levelshifts) {
169  levelshifts_ = levelshifts;
170  numLevels_ = levelshifts_.size();
171 }
172 
173 // initialize
174 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
175 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::initialize() {
176  TentPfact_ = rcp(new TentativePFactory);
177  Pfact_ = rcp(new SaPFactory);
178  PgPfact_ = rcp(new PgPFactory);
179  TransPfact_ = rcp(new TransPFactory);
180  Rfact_ = rcp(new GenericRFactory);
181  Acfact_ = rcp(new RAPFactory);
182  Acshift_ = rcp(new RAPShiftFactory);
183  Amalgfact_ = rcp(new AmalgamationFactory);
184  Dropfact_ = rcp(new CoalesceDropFactory);
185  UCaggfact_ = rcp(new UncoupledAggregationFactory);
186  CoarseMapfact_ = rcp(new CoarseMapFactory);
187  Manager_ = rcp(new FactoryManager);
188  Manager_->SetFactory("UnAmalgamationInfo", Amalgfact_);
189  Teuchos::ParameterList params;
190  params.set("lightweight wrap", true);
191  params.set("aggregation: drop scheme", "classical");
192  Dropfact_->SetParameterList(params);
193  Manager_->SetFactory("Graph", Dropfact_);
194  Manager_->SetFactory("Aggregates", UCaggfact_);
195  Manager_->SetFactory("CoarseMap", CoarseMapfact_);
196  Manager_->SetFactory("Ptent", TentPfact_);
197  if (isSymmetric_ == true) {
198  Manager_->SetFactory("P", Pfact_);
199  Manager_->SetFactory("R", TransPfact_);
200  } else {
201  Manager_->SetFactory("P", PgPfact_);
202  Manager_->SetFactory("R", Rfact_);
203  solverType_ = 10;
204  }
205 
206  // choose smoother
207  if (Smoother_ == "jacobi") {
208  precType_ = "RELAXATION";
209  precList_.set("relaxation: type", "Jacobi");
210  precList_.set("relaxation: sweeps", smoother_sweeps_);
211  precList_.set("relaxation: damping factor", smoother_damping_);
212  } else if (Smoother_ == "gauss-seidel") {
213  precType_ = "RELAXATION";
214  precList_.set("relaxation: type", "Gauss-Seidel");
215  precList_.set("relaxation: sweeps", smoother_sweeps_);
216  precList_.set("relaxation: damping factor", smoother_damping_);
217  } else if (Smoother_ == "symmetric gauss-seidel") {
218  precType_ = "RELAXATION";
219  precList_.set("relaxation: type", "Symmetric Gauss-Seidel");
220  precList_.set("relaxation: sweeps", smoother_sweeps_);
221  precList_.set("relaxation: damping factor", smoother_damping_);
222  } else if (Smoother_ == "chebyshev") {
223  precType_ = "CHEBYSHEV";
224  } else if (Smoother_ == "krylov") {
225  precType_ = "KRYLOV";
226  precList_.set("krylov: iteration type", krylov_type_);
227  precList_.set("krylov: number of iterations", krylov_iterations_);
228  precList_.set("krylov: residual tolerance", 1.0e-8);
229  precList_.set("krylov: block size", 1);
230  precList_.set("krylov: preconditioner type", krylov_preconditioner_);
231  precList_.set("relaxation: sweeps", 1);
232  solverType_ = 10;
233  } else if (Smoother_ == "ilut") {
234  precType_ = "ILUT";
235  precList_.set("fact: ilut level-of-fill", ilu_leveloffill_);
236  precList_.set("fact: absolute threshold", ilu_abs_thresh_);
237  precList_.set("fact: relative threshold", ilu_rel_thresh_);
238  precList_.set("fact: drop tolerance", ilu_drop_tol_);
239  precList_.set("fact: relax value", ilu_relax_val_);
240  } else if (Smoother_ == "riluk") {
241  precType_ = "RILUK";
242  precList_.set("fact: iluk level-of-fill", ilu_leveloffill_);
243  precList_.set("fact: absolute threshold", ilu_abs_thresh_);
244  precList_.set("fact: relative threshold", ilu_rel_thresh_);
245  precList_.set("fact: drop tolerance", ilu_drop_tol_);
246  precList_.set("fact: relax value", ilu_relax_val_);
247  } else if (Smoother_ == "schwarz") {
248  precType_ = "SCHWARZ";
249  precList_.set("schwarz: overlap level", schwarz_overlap_);
250  precList_.set("schwarz: combine mode", schwarz_combinemode_);
251  precList_.set("schwarz: use reordering", schwarz_usereorder_);
252  // precList_.set("schwarz: filter singletons", true); // Disabled due to issues w/ Ifpack2/Zoltan2 w.r.t. Issue #560 - CMS 8/26/16
253  precList_.set("order_method", schwarz_ordermethod_);
254  precList_.sublist("schwarz: reordering list").set("order_method", schwarz_ordermethod_);
255  precList_.sublist("schwarz: subdomain solver parameters").set("fact: ilut level-of-fill", ilu_leveloffill_);
256  precList_.sublist("schwarz: subdomain solver parameters").set("fact: absolute threshold", ilu_abs_thresh_);
257  precList_.sublist("schwarz: subdomain solver parameters").set("fact: relative threshold", ilu_rel_thresh_);
258  precList_.sublist("schwarz: subdomain solver parameters").set("fact: drop tolerance", ilu_drop_tol_);
259  precList_.sublist("schwarz: subdomain solver parameters").set("fact: relax value", ilu_relax_val_);
260  } else if (Smoother_ == "superilu") {
261  precType_ = "superlu";
262  precList_.set("RowPerm", ilu_rowperm_);
263  precList_.set("ColPerm", ilu_colperm_);
264  precList_.set("DiagPivotThresh", ilu_diagpivotthresh_);
265  precList_.set("ILU_DropRule", ilu_drop_rule_);
266  precList_.set("ILU_DropTol", ilu_drop_tol_);
267  precList_.set("ILU_FillFactor", ilu_leveloffill_);
268  precList_.set("ILU_Norm", ilu_normtype_);
269  precList_.set("ILU_MILU", ilu_milutype_);
270  precList_.set("ILU_FillTol", ilu_fill_tol_);
271  precList_.set("ILU_Flag", true);
272  } else if (Smoother_ == "superlu") {
273  precType_ = "superlu";
274  precList_.set("ColPerm", ilu_colperm_);
275  precList_.set("DiagPivotThresh", ilu_diagpivotthresh_);
276  }
277 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
278  // construct smoother
279  smooProto_ = rcp(new Ifpack2Smoother(precType_, precList_));
280  smooFact_ = rcp(new SmootherFactory(smooProto_));
281 #if defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLU)
282  coarsestSmooProto_ = rcp(new DirectSolver("Superlu", coarsestSmooList_));
283 #elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_KLU2)
284  coarsestSmooProto_ = rcp(new DirectSolver("Klu", coarsestSmooList_));
285 #elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLUDIST)
286  coarsestSmooProto_ = rcp(new DirectSolver("Superludist", coarsestSmooList_));
287 #else
288  coarsestSmooProto_ = rcp(new Ifpack2Smoother(precType_, precList_));
289 #endif
290  coarsestSmooFact_ = rcp(new SmootherFactory(coarsestSmooProto_, Teuchos::null));
291 
292  // For setupSlowRAP and setupFastRAP, the prolongation/restriction matrices
293  // are constructed with the stiffness matrix. These matrices are kept for future
294  // setup calls; this is achieved by calling Hierarchy->Keep(). It is particularly
295  // useful for multiple frequency problems - when the frequency/preconditioner
296  // changes, you only compute coarse grids (RAPs) and setup level smoothers when
297  // you call Hierarchy->Setup().
298  if (K_ != Teuchos::null) {
299  Manager_->SetFactory("Smoother", Teuchos::null);
300  Manager_->SetFactory("CoarseSolver", Teuchos::null);
301  Hierarchy_ = rcp(new Hierarchy(K_));
302  if (NullSpace_ != Teuchos::null)
303  Hierarchy_->GetLevel(0)->Set("Nullspace", NullSpace_);
304  if (isSymmetric_ == true) {
305  Hierarchy_->Keep("P", Pfact_.get());
306  Hierarchy_->Keep("R", TransPfact_.get());
307  Hierarchy_->SetImplicitTranspose(true);
308  } else {
309  Hierarchy_->Keep("P", PgPfact_.get());
310  Hierarchy_->Keep("R", Rfact_.get());
311  }
312  Hierarchy_->Keep("Ptent", TentPfact_.get());
313  Hierarchy_->SetMaxCoarseSize(coarseGridSize_);
314  Hierarchy_->Setup(*Manager_, 0, numLevels_);
315  GridTransfersExist_ = true;
316  }
317  // Use preconditioning matrix to setup prolongation/restriction operators
318  else {
319  Manager_->SetFactory("Smoother", smooFact_);
320  Manager_->SetFactory("CoarseSolver", coarsestSmooFact_);
321  Hierarchy_ = rcp(new Hierarchy(P_));
322  if (NullSpace_ != Teuchos::null)
323  Hierarchy_->GetLevel(0)->Set("Nullspace", NullSpace_);
324  if (isSymmetric_ == true)
325  Hierarchy_->SetImplicitTranspose(true);
326  Hierarchy_->SetMaxCoarseSize(coarseGridSize_);
327  Hierarchy_->Setup(*Manager_, 0, numLevels_);
328  GridTransfersExist_ = true;
329  }
330 
331  // Belos Linear Problem and Solver Manager
332  BelosList_ = rcp(new Teuchos::ParameterList("GMRES"));
333  BelosList_->set("Maximum Iterations", iters_);
334  BelosList_->set("Convergence Tolerance", tol_);
335  BelosList_->set("Verbosity", Belos::Errors + Belos::Warnings + Belos::StatusTestDetails);
336  BelosList_->set("Output Frequency", 1);
337  BelosList_->set("Output Style", Belos::Brief);
338  BelosList_->set("Num Blocks", restart_size_);
339  BelosList_->set("Num Recycled Blocks", recycle_size_);
340 #else
341  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "ShiftedLaplacian only available with Tpetra and GO=int enabled.");
342 #endif
343 }
344 
345 // setup coarse grids for new frequency
346 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
347 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupFastRAP() {
348  int numLevels = Hierarchy_->GetNumLevels();
349  Manager_->SetFactory("Smoother", smooFact_);
350  Manager_->SetFactory("CoarseSolver", coarsestSmooFact_);
351  Hierarchy_->GetLevel(0)->Set("A", P_);
352  Hierarchy_->Setup(*Manager_, 0, numLevels);
353  setupSolver();
354 }
355 
356 // setup coarse grids for new frequency
357 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
358 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupSlowRAP() {
359  int numLevels = Hierarchy_->GetNumLevels();
360  Acshift_->SetShifts(levelshifts_);
361  Manager_->SetFactory("Smoother", smooFact_);
362  Manager_->SetFactory("CoarseSolver", coarsestSmooFact_);
363  Manager_->SetFactory("A", Acshift_);
364  Manager_->SetFactory("K", Acshift_);
365  Manager_->SetFactory("M", Acshift_);
366  Hierarchy_->GetLevel(0)->Set("A", P_);
367  Hierarchy_->GetLevel(0)->Set("K", K_);
368  Hierarchy_->GetLevel(0)->Set("M", M_);
369  Hierarchy_->Setup(*Manager_, 0, numLevels);
370  setupSolver();
371 }
372 
373 // setup coarse grids for new frequency
374 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
375 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupNormalRAP() {
376  // Only setup hierarchy again if preconditioning matrix has changed
377  if (GridTransfersExist_ == false) {
378  Hierarchy_ = rcp(new Hierarchy(P_));
379  if (NullSpace_ != Teuchos::null)
380  Hierarchy_->GetLevel(0)->Set("Nullspace", NullSpace_);
381  if (isSymmetric_ == true)
382  Hierarchy_->SetImplicitTranspose(true);
383  Hierarchy_->SetMaxCoarseSize(coarseGridSize_);
384  Hierarchy_->Setup(*Manager_, 0, numLevels_);
385  GridTransfersExist_ = true;
386  }
387  setupSolver();
388 }
389 
390 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
391 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupSolver() {
392 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
393  // Define Preconditioner and Operator
394  MueLuOp_ = rcp(new MueLu::ShiftedLaplacianOperator<SC, LO, GO, NO>(Hierarchy_, A_, ncycles_, subiters_, option_, tol_));
395  // Belos Linear Problem
396  if (LinearProblem_ == Teuchos::null)
397  LinearProblem_ = rcp(new LinearProblem);
398  LinearProblem_->setOperator(TpetraA_);
399  LinearProblem_->setRightPrec(MueLuOp_);
400  if (SolverManager_ == Teuchos::null) {
401  std::string solverName;
402  SolverFactory_ = rcp(new SolverFactory());
403  if (solverType_ == 1) {
404  solverName = "Block GMRES";
405  } else if (solverType_ == 2) {
406  solverName = "Recycling GMRES";
407  } else {
408  solverName = "Flexible GMRES";
409  }
410  SolverManager_ = SolverFactory_->create(solverName, BelosList_);
411  SolverManager_->setProblem(LinearProblem_);
412  }
413 #else
414  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "ShiftedLaplacian only available with Tpetra and GO=int enabled.");
415 #endif
416 }
417 
418 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
419 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::resetLinearProblem() {
420 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
421  LinearProblem_->setOperator(TpetraA_);
422 #else
423  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "ShiftedLaplacian only available with Tpetra and GO=int enabled.");
424 #endif
425 }
426 
427 // Solve phase
428 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
429 int ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::solve(const RCP<TMV> B, RCP<TMV>& X) {
430 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
431  // Set left and right hand sides for Belos
432  LinearProblem_->setProblem(X, B);
433  // iterative solve
434  SolverManager_->solve();
435 #else
436  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "ShiftedLaplacian only available with Tpetra and GO=int enabled.");
437 #endif
438  return 0;
439 }
440 
441 // Solve phase
442 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
443 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multigrid_apply(const RCP<MultiVector> B,
444  RCP<MultiVector>& X) {
445  // Set left and right hand sides for Belos
446  Hierarchy_->Iterate(*B, *X, 1, true, 0);
447 }
448 
449 // Solve phase
450 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
451 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multigrid_apply(const RCP<Tpetra::MultiVector<SC, LO, GO, NO> > B,
452  RCP<Tpetra::MultiVector<SC, LO, GO, NO> >& X) {
453  Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > XpetraX = Teuchos::rcp(new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(X));
454  Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > XpetraB = Teuchos::rcp(new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(B));
455  // Set left and right hand sides for Belos
456  Hierarchy_->Iterate(*XpetraB, *XpetraX, 1, true, 0);
457 }
458 
459 // Get most recent iteration count
460 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
461 int ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetIterations() {
462 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
463  int numiters = SolverManager_->getNumIters();
464  return numiters;
465 #else
466  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "ShiftedLaplacian only available with Tpetra and GO=int enabled.");
467  return -1;
468 #endif
469 }
470 
471 // Get most recent solver tolerance achieved
472 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
473 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
474 ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetResidual() {
475  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
476 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
477  MT residual = SolverManager_->achievedTol();
478  return residual;
479 #else
480  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "ShiftedLaplacian only available with Tpetra and GO=int enabled.");
481  return MT(-1.0);
482 #endif
483 }
484 
485 } // namespace MueLu
486 
487 #define MUELU_SHIFTEDLAPLACIAN_SHORT
488 
489 #endif // if defined(HAVE_MUELU_IFPACK2)
490 #endif // MUELU_SHIFTEDLAPLACIAN_DEF_HPP
MueLu::SmootherFactory
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
Definition: MueLu_SmootherFactory_decl.hpp:54
MueLu::FactoryManager
This class specifies the default factory that should generate some data on a Level if the data does n...
Definition: MueLu_FactoryManager_decl.hpp:70
MueLu::CoarseMapFactory
Factory for generating coarse level map. Used by TentativePFactory.
Definition: MueLu_CoarseMapFactory_decl.hpp:77
MueLu::RAPShiftFactory
Factory for building coarse grid matrices, when the matrix is of the form K+a*M. Useful when you want...
Definition: MueLu_RAPShiftFactory_decl.hpp:40
MueLu::ShiftedLaplacian::multigrid_apply
void multigrid_apply(const RCP< MultiVector > B, RCP< MultiVector > &X)
Definition: MueLu_ShiftedLaplacian_def.hpp:443
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
MueLu::TentativePFactory
Factory for building tentative prolongator.
Definition: MueLu_TentativePFactory_decl.hpp:85
MueLu::ShiftedLaplacian::setParameters
void setParameters(Teuchos::RCP< Teuchos::ParameterList > paramList)
Definition: MueLu_ShiftedLaplacian_def.hpp:47
MueLu::DirectSolver
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
Definition: MueLu_DirectSolver_decl.hpp:43
MueLu::GenericRFactory
Factory for building restriction operators using a prolongator factory.
Definition: MueLu_GenericRFactory_decl.hpp:38
MueLu::ShiftedLaplacian::GetResidual
Teuchos::ScalarTraits< Scalar >::magnitudeType GetResidual()
Definition: MueLu_ShiftedLaplacian_def.hpp:474
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::ShiftedLaplacian::setLevelShifts
void setLevelShifts(std::vector< Scalar > levelshifts)
Definition: MueLu_ShiftedLaplacian_def.hpp:168
MueLu::AmalgamationFactory
AmalgamationFactory for subblocks of strided map based amalgamation data.
Definition: MueLu_AmalgamationFactory_decl.hpp:42
MueLu_ShiftedLaplacian_decl.hpp
MueLu::ShiftedLaplacian::setcoords
void setcoords(RCP< MultiVector > &Coords)
Definition: MueLu_ShiftedLaplacian_def.hpp:158
MueLu::ShiftedLaplacianOperator
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator, with an optional two-level correction...
Definition: MueLu_ShiftedLaplacianOperator_decl.hpp:32
MueLu::ShiftedLaplacian::solve
int solve(const RCP< TMV > B, RCP< TMV > &X)
Definition: MueLu_ShiftedLaplacian_def.hpp:429
MueLu::CoalesceDropFactory
Factory for creating a graph based on a given matrix.
Definition: MueLu_CoalesceDropFactory_decl.hpp:98
MueLu::ShiftedLaplacian::setProblemMatrix
void setProblemMatrix(RCP< Matrix > &A)
Definition: MueLu_ShiftedLaplacian_def.hpp:101
MueLu::TransPFactory
Factory for building restriction operators.
Definition: MueLu_TransPFactory_decl.hpp:37
Teuchos::ParameterList::sublist
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
residual
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:36
MueLu::Warnings
Print all warning messages.
Definition: MueLu_VerbosityLevel.hpp:52
MueLu::RAPFactory
Factory for building coarse matrices.
Definition: MueLu_RAPFactory_decl.hpp:38
MueLu::PgPFactory
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
Definition: MueLu_PgPFactory_decl.hpp:46
MueLu::Ifpack2Smoother
Class that encapsulates Ifpack2 smoothers.
Definition: MueLu_Ifpack2Smoother_decl.hpp:55
MueLu::SaPFactory
Factory for building Smoothed Aggregation prolongators.Input/output of SaPFactory
Definition: MueLu_SaPFactory_decl.hpp:62
MueLu::UncoupledAggregationFactory
Factory for building uncoupled aggregates.
Definition: MueLu_UncoupledAggregationFactory_decl.hpp:104
MueLu::Utilities::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
Definition: MueLu_Utilities_def.hpp:224
MueLu::ShiftedLaplacian::setNullSpace
void setNullSpace(RCP< MultiVector > NullSpace)
Definition: MueLu_ShiftedLaplacian_def.hpp:163
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition: MueLu_Hierarchy_decl.hpp:63
MueLu::ShiftedLaplacian::setPreconditioningMatrix
void setPreconditioningMatrix(RCP< Matrix > &P)
Definition: MueLu_ShiftedLaplacian_def.hpp:123