MueLu Version of the Day
Loading...
Searching...
No Matches
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
39namespace MueLu {
40
41// Destructor
42template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
43ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~ShiftedLaplacian() {}
44
45// Input
46template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
47void 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
100template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
101void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setProblemMatrix(RCP<Matrix>& A) {
102 A_ = A;
103 if (A_ != Teuchos::null)
104 TpetraA_ = toTpetra(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
113template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
114void 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
122template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
123void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setPreconditioningMatrix(RCP<Matrix>& P) {
124 P_ = P;
125 GridTransfersExist_ = false;
126}
127
128template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
129void 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
135template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
139
140template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141void 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
146template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
150
151template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
152void 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
157template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
158void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setcoords(RCP<MultiVector>& Coords) {
159 Coords_ = Coords;
160}
161
162template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
163void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setNullSpace(RCP<MultiVector> NullSpace) {
164 NullSpace_ = NullSpace;
165}
166
167template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
168void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setLevelShifts(std::vector<Scalar> levelshifts) {
169 levelshifts_ = levelshifts;
170 numLevels_ = levelshifts_.size();
171}
172
173// initialize
174template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
175void 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
346template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
347void 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
357template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
358void 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
374template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
375void 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
390template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
391void 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
418template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
419void 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
428template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
429int 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
442template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
443void 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
450template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
451void 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
460template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
461int 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
472template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
473typename Teuchos::ScalarTraits<Scalar>::magnitudeType
474ShiftedLaplacian<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_ShiftedLaplacian_decl.hpp
MueLu::AmalgamationFactory
AmalgamationFactory for subblocks of strided map based amalgamation data.
Definition MueLu_AmalgamationFactory_decl.hpp:42
MueLu::CoalesceDropFactory
Factory for creating a graph based on a given matrix.
Definition MueLu_CoalesceDropFactory_decl.hpp:98
MueLu::CoarseMapFactory
Factory for generating coarse level map. Used by TentativePFactory.
Definition MueLu_CoarseMapFactory_decl.hpp:77
MueLu::DirectSolver
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
Definition MueLu_DirectSolver_decl.hpp:43
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition MueLu_Exceptions.hpp:36
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::GenericRFactory
Factory for building restriction operators using a prolongator factory.
Definition MueLu_GenericRFactory_decl.hpp:38
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition MueLu_Hierarchy_decl.hpp:63
MueLu::Ifpack2Smoother
Class that encapsulates Ifpack2 smoothers.
Definition MueLu_Ifpack2Smoother_decl.hpp:55
MueLu::PgPFactory
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
Definition MueLu_PgPFactory_decl.hpp:46
MueLu::RAPFactory
Factory for building coarse matrices.
Definition MueLu_RAPFactory_decl.hpp:38
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::SaPFactory
Factory for building Smoothed Aggregation prolongators.
Definition MueLu_SaPFactory_decl.hpp:62
MueLu::ShiftedLaplacianOperator
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator, with an optional two-level correction....
Definition MueLu_ShiftedLaplacianOperator_decl.hpp:33
MueLu::ShiftedLaplacian::setPreconditioningMatrix
void setPreconditioningMatrix(RCP< Matrix > &P)
Definition MueLu_ShiftedLaplacian_def.hpp:123
MueLu::ShiftedLaplacian::setNullSpace
void setNullSpace(RCP< MultiVector > NullSpace)
Definition MueLu_ShiftedLaplacian_def.hpp:163
MueLu::ShiftedLaplacian::setcoords
void setcoords(RCP< MultiVector > &Coords)
Definition MueLu_ShiftedLaplacian_def.hpp:158
MueLu::ShiftedLaplacian::GetResidual
Teuchos::ScalarTraits< Scalar >::magnitudeType GetResidual()
Definition MueLu_ShiftedLaplacian_def.hpp:474
MueLu::ShiftedLaplacian::solve
int solve(const RCP< TMV > B, RCP< TMV > &X)
Definition MueLu_ShiftedLaplacian_def.hpp:429
MueLu::ShiftedLaplacian::setLevelShifts
void setLevelShifts(std::vector< Scalar > levelshifts)
Definition MueLu_ShiftedLaplacian_def.hpp:168
MueLu::ShiftedLaplacian::setParameters
void setParameters(Teuchos::RCP< Teuchos::ParameterList > paramList)
Definition MueLu_ShiftedLaplacian_def.hpp:47
MueLu::ShiftedLaplacian::setProblemMatrix
void setProblemMatrix(RCP< Matrix > &A)
Definition MueLu_ShiftedLaplacian_def.hpp:101
MueLu::ShiftedLaplacian::multigrid_apply
void multigrid_apply(const RCP< MultiVector > B, RCP< MultiVector > &X)
Definition MueLu_ShiftedLaplacian_def.hpp:443
MueLu::SmootherFactory
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
Definition MueLu_SmootherFactory_decl.hpp:54
MueLu::TentativePFactory
Factory for building tentative prolongator.
Definition MueLu_TentativePFactory_decl.hpp:85
MueLu::TransPFactory
Factory for building restriction operators.
Definition MueLu_TransPFactory_decl.hpp:37
MueLu
Namespace for MueLu classes and methods.
Definition MueLu_BrickAggregationFactory_decl.hpp:42