MueLu Version of the Day
Loading...
Searching...
No Matches
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_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
76namespace MueLu {
77
78// Destructor
79template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
81
82// Input
83template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84void 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
124template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
139template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
140void 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
150template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
152
153 P_=P;
154 GridTransfersExist_=false;
155
156}
157
158template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
168template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
170
171 K_=K;
172
173}
174
175template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
176void 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
184template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
186
187 M_=M;
188
189}
190
191template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
192void 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
200template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
202
203 Coords_=Coords;
204
205}
206
207template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
209
210 NullSpace_=NullSpace;
211
212}
213
214template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
216
217 levelshifts_=levelshifts;
218 numLevels_=levelshifts_.size();
219
220}
221
222// initialize
223template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
413template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
426template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
445template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
463template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
465
466#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
467 // Define Preconditioner and Operator
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
489template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
500template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
515template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
524template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
525void 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
537template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
550template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
551typename Teuchos::ScalarTraits<Scalar>::magnitudeType
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
AmalgamationFactory for subblocks of strided map based amalgamation data.
Factory for creating a graph based on a given matrix.
Factory for generating coarse level map. Used by TentativePFactory.
Factory for coarsening a graph with uncoupled aggregation.
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
Exception throws to report errors in the internal logical of the program.
This class specifies the default factory that should generate some data on a Level if the data does n...
Factory for building restriction operators using a prolongator factory.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Class that encapsulates Ifpack2 smoothers.
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
Factory for building coarse matrices.
Factory for building coarse grid matrices, when the matrix is of the form K+a*M. Useful when you want...
Factory for building Smoothed Aggregation prolongators.
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator, with an optional two-level correction....
void setPreconditioningMatrix(RCP< Matrix > &P)
void setNullSpace(RCP< MultiVector > NullSpace)
void setcoords(RCP< MultiVector > &Coords)
Teuchos::ScalarTraits< Scalar >::magnitudeType GetResidual()
int solve(const RCP< TMV > B, RCP< TMV > &X)
void setLevelShifts(std::vector< Scalar > levelshifts)
void setParameters(Teuchos::RCP< Teuchos::ParameterList > paramList)
void setProblemMatrix(RCP< Matrix > &A)
void multigrid_apply(const RCP< MultiVector > B, RCP< MultiVector > &X)
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
Factory for building tentative prolongator.
Factory for building restriction operators.
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Namespace for MueLu classes and methods.
@ Keep
Always keep data, even accross run. This flag is set by Level::Keep(). This flag is propagated to coa...