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