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