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