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