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