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_ = toTpetra(A_);
105  if (LinearProblem_ != Teuchos::null)
106  LinearProblem_->setOperator(TpetraA_);
107 }
108 
109 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
110 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setProblemMatrix(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraA) {
111  TpetraA_ = TpetraA;
112  if (LinearProblem_ != Teuchos::null)
113  LinearProblem_->setOperator(TpetraA_);
114 }
115 
116 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
117 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setPreconditioningMatrix(RCP<Matrix>& P) {
118  P_ = P;
119  GridTransfersExist_ = false;
120 }
121 
122 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
123 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setPreconditioningMatrix(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraP) {
124  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Atmp = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraP));
125  P_ = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp));
126  GridTransfersExist_ = false;
127 }
128 
129 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
130 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setstiff(RCP<Matrix>& K) {
131  K_ = K;
132 }
133 
134 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
135 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setstiff(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraK) {
136  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Atmp = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraK));
137  K_ = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp));
138 }
139 
140 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setmass(RCP<Matrix>& M) {
142  M_ = M;
143 }
144 
145 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
146 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setmass(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraM) {
147  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Atmp = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraM));
148  M_ = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp));
149 }
150 
151 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
152 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setcoords(RCP<MultiVector>& Coords) {
153  Coords_ = Coords;
154 }
155 
156 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
157 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setNullSpace(RCP<MultiVector> NullSpace) {
158  NullSpace_ = NullSpace;
159 }
160 
161 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
162 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setLevelShifts(std::vector<Scalar> levelshifts) {
163  levelshifts_ = levelshifts;
164  numLevels_ = levelshifts_.size();
165 }
166 
167 // initialize
168 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
169 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::initialize() {
170  TentPfact_ = rcp(new TentativePFactory);
171  Pfact_ = rcp(new SaPFactory);
172  PgPfact_ = rcp(new PgPFactory);
173  TransPfact_ = rcp(new TransPFactory);
174  Rfact_ = rcp(new GenericRFactory);
175  Acfact_ = rcp(new RAPFactory);
176  Acshift_ = rcp(new RAPShiftFactory);
177  Amalgfact_ = rcp(new AmalgamationFactory);
178  Dropfact_ = rcp(new CoalesceDropFactory);
179  UCaggfact_ = rcp(new UncoupledAggregationFactory);
180  CoarseMapfact_ = rcp(new CoarseMapFactory);
181  Manager_ = rcp(new FactoryManager);
182  Manager_->SetFactory("UnAmalgamationInfo", Amalgfact_);
183  Teuchos::ParameterList params;
184  params.set("lightweight wrap", true);
185  params.set("aggregation: drop scheme", "classical");
186  Dropfact_->SetParameterList(params);
187  Manager_->SetFactory("Graph", Dropfact_);
188  Manager_->SetFactory("Aggregates", UCaggfact_);
189  Manager_->SetFactory("CoarseMap", CoarseMapfact_);
190  Manager_->SetFactory("Ptent", TentPfact_);
191  if (isSymmetric_ == true) {
192  Manager_->SetFactory("P", Pfact_);
193  Manager_->SetFactory("R", TransPfact_);
194  } else {
195  Manager_->SetFactory("P", PgPfact_);
196  Manager_->SetFactory("R", Rfact_);
197  solverType_ = 10;
198  }
199 
200  // choose smoother
201  if (Smoother_ == "jacobi") {
202  precType_ = "RELAXATION";
203  precList_.set("relaxation: type", "Jacobi");
204  precList_.set("relaxation: sweeps", smoother_sweeps_);
205  precList_.set("relaxation: damping factor", smoother_damping_);
206  } else if (Smoother_ == "gauss-seidel") {
207  precType_ = "RELAXATION";
208  precList_.set("relaxation: type", "Gauss-Seidel");
209  precList_.set("relaxation: sweeps", smoother_sweeps_);
210  precList_.set("relaxation: damping factor", smoother_damping_);
211  } else if (Smoother_ == "symmetric gauss-seidel") {
212  precType_ = "RELAXATION";
213  precList_.set("relaxation: type", "Symmetric Gauss-Seidel");
214  precList_.set("relaxation: sweeps", smoother_sweeps_);
215  precList_.set("relaxation: damping factor", smoother_damping_);
216  } else if (Smoother_ == "chebyshev") {
217  precType_ = "CHEBYSHEV";
218  } else if (Smoother_ == "krylov") {
219  precType_ = "KRYLOV";
220  precList_.set("krylov: iteration type", krylov_type_);
221  precList_.set("krylov: number of iterations", krylov_iterations_);
222  precList_.set("krylov: residual tolerance", 1.0e-8);
223  precList_.set("krylov: block size", 1);
224  precList_.set("krylov: preconditioner type", krylov_preconditioner_);
225  precList_.set("relaxation: sweeps", 1);
226  solverType_ = 10;
227  } else if (Smoother_ == "ilut") {
228  precType_ = "ILUT";
229  precList_.set("fact: ilut level-of-fill", ilu_leveloffill_);
230  precList_.set("fact: absolute threshold", ilu_abs_thresh_);
231  precList_.set("fact: relative threshold", ilu_rel_thresh_);
232  precList_.set("fact: drop tolerance", ilu_drop_tol_);
233  precList_.set("fact: relax value", ilu_relax_val_);
234  } else if (Smoother_ == "riluk") {
235  precType_ = "RILUK";
236  precList_.set("fact: iluk level-of-fill", ilu_leveloffill_);
237  precList_.set("fact: absolute threshold", ilu_abs_thresh_);
238  precList_.set("fact: relative threshold", ilu_rel_thresh_);
239  precList_.set("fact: drop tolerance", ilu_drop_tol_);
240  precList_.set("fact: relax value", ilu_relax_val_);
241  } else if (Smoother_ == "schwarz") {
242  precType_ = "SCHWARZ";
243  precList_.set("schwarz: overlap level", schwarz_overlap_);
244  precList_.set("schwarz: combine mode", schwarz_combinemode_);
245  precList_.set("schwarz: use reordering", schwarz_usereorder_);
246  // precList_.set("schwarz: filter singletons", true); // Disabled due to issues w/ Ifpack2/Zoltan2 w.r.t. Issue #560 - CMS 8/26/16
247  precList_.set("order_method", schwarz_ordermethod_);
248  precList_.sublist("schwarz: reordering list").set("order_method", schwarz_ordermethod_);
249  precList_.sublist("schwarz: subdomain solver parameters").set("fact: ilut level-of-fill", ilu_leveloffill_);
250  precList_.sublist("schwarz: subdomain solver parameters").set("fact: absolute threshold", ilu_abs_thresh_);
251  precList_.sublist("schwarz: subdomain solver parameters").set("fact: relative threshold", ilu_rel_thresh_);
252  precList_.sublist("schwarz: subdomain solver parameters").set("fact: drop tolerance", ilu_drop_tol_);
253  precList_.sublist("schwarz: subdomain solver parameters").set("fact: relax value", ilu_relax_val_);
254  } else if (Smoother_ == "superilu") {
255  precType_ = "superlu";
256  precList_.set("RowPerm", ilu_rowperm_);
257  precList_.set("ColPerm", ilu_colperm_);
258  precList_.set("DiagPivotThresh", ilu_diagpivotthresh_);
259  precList_.set("ILU_DropRule", ilu_drop_rule_);
260  precList_.set("ILU_DropTol", ilu_drop_tol_);
261  precList_.set("ILU_FillFactor", ilu_leveloffill_);
262  precList_.set("ILU_Norm", ilu_normtype_);
263  precList_.set("ILU_MILU", ilu_milutype_);
264  precList_.set("ILU_FillTol", ilu_fill_tol_);
265  precList_.set("ILU_Flag", true);
266  } else if (Smoother_ == "superlu") {
267  precType_ = "superlu";
268  precList_.set("ColPerm", ilu_colperm_);
269  precList_.set("DiagPivotThresh", ilu_diagpivotthresh_);
270  }
271  // construct smoother
272  smooProto_ = rcp(new Ifpack2Smoother(precType_, precList_));
273  smooFact_ = rcp(new SmootherFactory(smooProto_));
274 #if defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLU)
275  coarsestSmooProto_ = rcp(new DirectSolver("Superlu", coarsestSmooList_));
276 #elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_KLU2)
277  coarsestSmooProto_ = rcp(new DirectSolver("Klu", coarsestSmooList_));
278 #elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLUDIST)
279  coarsestSmooProto_ = rcp(new DirectSolver("Superludist", coarsestSmooList_));
280 #else
281  coarsestSmooProto_ = rcp(new Ifpack2Smoother(precType_, precList_));
282 #endif
283  coarsestSmooFact_ = rcp(new SmootherFactory(coarsestSmooProto_, Teuchos::null));
284 
285  // For setupSlowRAP and setupFastRAP, the prolongation/restriction matrices
286  // are constructed with the stiffness matrix. These matrices are kept for future
287  // setup calls; this is achieved by calling Hierarchy->Keep(). It is particularly
288  // useful for multiple frequency problems - when the frequency/preconditioner
289  // changes, you only compute coarse grids (RAPs) and setup level smoothers when
290  // you call Hierarchy->Setup().
291  if (K_ != Teuchos::null) {
292  Manager_->SetFactory("Smoother", Teuchos::null);
293  Manager_->SetFactory("CoarseSolver", Teuchos::null);
294  Hierarchy_ = rcp(new Hierarchy(K_));
295  if (NullSpace_ != Teuchos::null)
296  Hierarchy_->GetLevel(0)->Set("Nullspace", NullSpace_);
297  if (isSymmetric_ == true) {
298  Hierarchy_->Keep("P", Pfact_.get());
299  Hierarchy_->Keep("R", TransPfact_.get());
300  Hierarchy_->SetImplicitTranspose(true);
301  } else {
302  Hierarchy_->Keep("P", PgPfact_.get());
303  Hierarchy_->Keep("R", Rfact_.get());
304  }
305  Hierarchy_->Keep("Ptent", TentPfact_.get());
306  Hierarchy_->SetMaxCoarseSize(coarseGridSize_);
307  Hierarchy_->Setup(*Manager_, 0, numLevels_);
308  GridTransfersExist_ = true;
309  }
310  // Use preconditioning matrix to setup prolongation/restriction operators
311  else {
312  Manager_->SetFactory("Smoother", smooFact_);
313  Manager_->SetFactory("CoarseSolver", coarsestSmooFact_);
314  Hierarchy_ = rcp(new Hierarchy(P_));
315  if (NullSpace_ != Teuchos::null)
316  Hierarchy_->GetLevel(0)->Set("Nullspace", NullSpace_);
317  if (isSymmetric_ == true)
318  Hierarchy_->SetImplicitTranspose(true);
319  Hierarchy_->SetMaxCoarseSize(coarseGridSize_);
320  Hierarchy_->Setup(*Manager_, 0, numLevels_);
321  GridTransfersExist_ = true;
322  }
323 
324  // Belos Linear Problem and Solver Manager
325  BelosList_ = rcp(new Teuchos::ParameterList("GMRES"));
326  BelosList_->set("Maximum Iterations", iters_);
327  BelosList_->set("Convergence Tolerance", tol_);
328  BelosList_->set("Verbosity", Belos::Errors + Belos::Warnings + Belos::StatusTestDetails);
329  BelosList_->set("Output Frequency", 1);
330  BelosList_->set("Output Style", Belos::Brief);
331  BelosList_->set("Num Blocks", restart_size_);
332  BelosList_->set("Num Recycled Blocks", recycle_size_);
333 }
334 
335 // setup coarse grids for new frequency
336 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
337 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupFastRAP() {
338  int numLevels = Hierarchy_->GetNumLevels();
339  Manager_->SetFactory("Smoother", smooFact_);
340  Manager_->SetFactory("CoarseSolver", coarsestSmooFact_);
341  Hierarchy_->GetLevel(0)->Set("A", P_);
342  Hierarchy_->Setup(*Manager_, 0, numLevels);
343  setupSolver();
344 }
345 
346 // setup coarse grids for new frequency
347 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
348 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupSlowRAP() {
349  int numLevels = Hierarchy_->GetNumLevels();
350  Acshift_->SetShifts(levelshifts_);
351  Manager_->SetFactory("Smoother", smooFact_);
352  Manager_->SetFactory("CoarseSolver", coarsestSmooFact_);
353  Manager_->SetFactory("A", Acshift_);
354  Manager_->SetFactory("K", Acshift_);
355  Manager_->SetFactory("M", Acshift_);
356  Hierarchy_->GetLevel(0)->Set("A", P_);
357  Hierarchy_->GetLevel(0)->Set("K", K_);
358  Hierarchy_->GetLevel(0)->Set("M", M_);
359  Hierarchy_->Setup(*Manager_, 0, numLevels);
360  setupSolver();
361 }
362 
363 // setup coarse grids for new frequency
364 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
365 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupNormalRAP() {
366  // Only setup hierarchy again if preconditioning matrix has changed
367  if (GridTransfersExist_ == false) {
368  Hierarchy_ = rcp(new Hierarchy(P_));
369  if (NullSpace_ != Teuchos::null)
370  Hierarchy_->GetLevel(0)->Set("Nullspace", NullSpace_);
371  if (isSymmetric_ == true)
372  Hierarchy_->SetImplicitTranspose(true);
373  Hierarchy_->SetMaxCoarseSize(coarseGridSize_);
374  Hierarchy_->Setup(*Manager_, 0, numLevels_);
375  GridTransfersExist_ = true;
376  }
377  setupSolver();
378 }
379 
380 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
381 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupSolver() {
382  // Define Preconditioner and Operator
383  MueLuOp_ = rcp(new MueLu::ShiftedLaplacianOperator<SC, LO, GO, NO>(Hierarchy_, A_, ncycles_, subiters_, option_, tol_));
384  // Belos Linear Problem
385  if (LinearProblem_ == Teuchos::null)
386  LinearProblem_ = rcp(new LinearProblem);
387  LinearProblem_->setOperator(TpetraA_);
388  LinearProblem_->setRightPrec(MueLuOp_);
389  if (SolverManager_ == Teuchos::null) {
390  std::string solverName;
391  SolverFactory_ = rcp(new SolverFactory());
392  if (solverType_ == 1) {
393  solverName = "Block GMRES";
394  } else if (solverType_ == 2) {
395  solverName = "Recycling GMRES";
396  } else {
397  solverName = "Flexible GMRES";
398  }
399  SolverManager_ = SolverFactory_->create(solverName, BelosList_);
400  SolverManager_->setProblem(LinearProblem_);
401  }
402 }
403 
404 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
405 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::resetLinearProblem() {
406  LinearProblem_->setOperator(TpetraA_);
407 }
408 
409 // Solve phase
410 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
411 Belos::ReturnType ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::solve(const RCP<TMV> B, RCP<TMV>& X) {
412  // Set left and right hand sides for Belos
413  LinearProblem_->setProblem(X, B);
414  // iterative solve
415  return SolverManager_->solve();
416 }
417 
418 // Solve phase
419 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
420 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multigrid_apply(const RCP<MultiVector> B,
421  RCP<MultiVector>& X) {
422  // Set left and right hand sides for Belos
423  Hierarchy_->Iterate(*B, *X, 1, true, 0);
424 }
425 
426 // Solve phase
427 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
428 void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multigrid_apply(const RCP<Tpetra::MultiVector<SC, LO, GO, NO> > B,
429  RCP<Tpetra::MultiVector<SC, LO, GO, NO> >& X) {
430  Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > XpetraX = Teuchos::rcp(new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(X));
431  Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > XpetraB = Teuchos::rcp(new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(B));
432  // Set left and right hand sides for Belos
433  Hierarchy_->Iterate(*XpetraB, *XpetraX, 1, true, 0);
434 }
435 
436 // Get most recent iteration count
437 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
438 int ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetIterations() {
439  int numiters = SolverManager_->getNumIters();
440  return numiters;
441 }
442 
443 // Get most recent solver tolerance achieved
444 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
445 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
446 ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetResidual() {
447  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
448  MT residual = SolverManager_->achievedTol();
449  return residual;
450 }
451 
452 } // namespace MueLu
453 
454 #define MUELU_SHIFTEDLAPLACIAN_SHORT
455 
456 #endif // if defined(HAVE_MUELU_IFPACK2)
457 #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:420
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
MueLu::ShiftedLaplacian::solve
Belos::ReturnType solve(const RCP< TMV > B, RCP< TMV > &X)
Definition: MueLu_ShiftedLaplacian_def.hpp:411
MueLu::ShiftedLaplacian::SolverFactory
Belos::SolverFactory< SC, TMV, OP > SolverFactory
Definition: MueLu_ShiftedLaplacian_decl.hpp:80
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:446
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:162
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:152
MueLu::ShiftedLaplacianOperator
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator, with an optional two-level correction...
Definition: MueLu_ShiftedLaplacianOperator_decl.hpp:32
toTpetra
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
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::Warnings
Print all warning messages.
Definition: MueLu_VerbosityLevel.hpp:52
MueLu::ShiftedLaplacian::LinearProblem
Belos::LinearProblem< SC, TMV, OP > LinearProblem
Definition: MueLu_ShiftedLaplacian_decl.hpp:78
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::ShiftedLaplacian::setNullSpace
void setNullSpace(RCP< MultiVector > NullSpace)
Definition: MueLu_ShiftedLaplacian_def.hpp:157
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:117