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_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>
80ShiftedLaplacian<Scalar,LocalOrdinal,GlobalOrdinal,Node>::~ShiftedLaplacian() {}
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>
125void 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
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>
157
158template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
159void 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
168template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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>
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>
201void ShiftedLaplacian<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setcoords(RCP<MultiVector>& Coords) {
202
203 Coords_=Coords;
204
205}
206
207template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
208void ShiftedLaplacian<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setNullSpace(RCP<MultiVector> NullSpace) {
209
210 NullSpace_=NullSpace;
211
212}
213
214template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
215void ShiftedLaplacian<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setLevelShifts(std::vector<Scalar> levelshifts) {
216
217 levelshifts_=levelshifts;
218 numLevels_=levelshifts_.size();
219
220}
221
222// initialize
223template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
224void 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
413template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
414void 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
426template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
427void 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
445template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
446void 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
463template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
464void 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
489template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
490void 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
500template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
501int 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
515template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
516void 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
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>
538int 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
550template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
551typename Teuchos::ScalarTraits<Scalar>::magnitudeType
552ShiftedLaplacian<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_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:136
MueLu::CoarseMapFactory
Factory for generating coarse level map. Used by TentativePFactory.
Definition MueLu_CoarseMapFactory_decl.hpp:114
MueLu::CoupledAggregationFactory
Factory for coarsening a graph with uncoupled aggregation.
Definition MueLu_CoupledAggregationFactory_decl.hpp:84
MueLu::DirectSolver
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
Definition MueLu_DirectSolver_decl.hpp:81
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:112
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:103
MueLu::Ifpack2Smoother
Class that encapsulates Ifpack2 smoothers.
Definition MueLu_Ifpack2Smoother_decl.hpp:94
MueLu::PgPFactory
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
Definition MueLu_PgPFactory_decl.hpp:84
MueLu::RAPFactory
Factory for building coarse matrices.
Definition MueLu_RAPFactory_decl.hpp:77
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::SaPFactory
Factory for building Smoothed Aggregation prolongators.
Definition MueLu_SaPFactory_decl.hpp:99
MueLu::ShiftedLaplacianOperator
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator, with an optional two-level correction....
Definition MueLu_ShiftedLaplacianOperator_decl.hpp:71
MueLu::ShiftedLaplacian::setPreconditioningMatrix
void setPreconditioningMatrix(RCP< Matrix > &P)
Definition MueLu_ShiftedLaplacian_def.hpp:151
MueLu::ShiftedLaplacian::setNullSpace
void setNullSpace(RCP< MultiVector > NullSpace)
Definition MueLu_ShiftedLaplacian_def.hpp:208
MueLu::ShiftedLaplacian::MueLuOp_
RCP< MueLu::ShiftedLaplacianOperator< SC, LO, GO, NO > > MueLuOp_
Definition MueLu_ShiftedLaplacian_decl.hpp:299
MueLu::ShiftedLaplacian::coarsestSmooProto_
RCP< SmootherPrototype > coarsestSmooProto_
Definition MueLu_ShiftedLaplacian_decl.hpp:292
MueLu::ShiftedLaplacian::setcoords
void setcoords(RCP< MultiVector > &Coords)
Definition MueLu_ShiftedLaplacian_def.hpp:201
MueLu::ShiftedLaplacian::GetResidual
Teuchos::ScalarTraits< Scalar >::magnitudeType GetResidual()
Definition MueLu_ShiftedLaplacian_def.hpp:552
MueLu::ShiftedLaplacian::TpetraA_
RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > TpetraA_
Definition MueLu_ShiftedLaplacian_decl.hpp:300
MueLu::ShiftedLaplacian::UCaggfact_
RCP< UncoupledAggregationFactory > UCaggfact_
Definition MueLu_ShiftedLaplacian_decl.hpp:290
MueLu::ShiftedLaplacian::solve
int solve(const RCP< TMV > B, RCP< TMV > &X)
Definition MueLu_ShiftedLaplacian_def.hpp:501
MueLu::ShiftedLaplacian::setLevelShifts
void setLevelShifts(std::vector< Scalar > levelshifts)
Definition MueLu_ShiftedLaplacian_def.hpp:215
MueLu::ShiftedLaplacian::setParameters
void setParameters(Teuchos::RCP< Teuchos::ParameterList > paramList)
Definition MueLu_ShiftedLaplacian_def.hpp:84
MueLu::ShiftedLaplacian::Aggfact_
RCP< CoupledAggregationFactory > Aggfact_
Definition MueLu_ShiftedLaplacian_decl.hpp:289
MueLu::ShiftedLaplacian::setProblemMatrix
void setProblemMatrix(RCP< Matrix > &A)
Definition MueLu_ShiftedLaplacian_def.hpp:125
MueLu::ShiftedLaplacian::multigrid_apply
void multigrid_apply(const RCP< MultiVector > B, RCP< MultiVector > &X)
Definition MueLu_ShiftedLaplacian_def.hpp:516
MueLu::ShiftedLaplacian::Amalgfact_
RCP< AmalgamationFactory > Amalgfact_
Definition MueLu_ShiftedLaplacian_decl.hpp:287
MueLu::ShiftedLaplacian::Dropfact_
RCP< CoalesceDropFactory > Dropfact_
Definition MueLu_ShiftedLaplacian_decl.hpp:288
MueLu::SmootherFactory
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
Definition MueLu_SmootherFactory_decl.hpp:91
MueLu::TentativePFactory
Factory for building tentative prolongator.
Definition MueLu_TentativePFactory_decl.hpp:123
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:263
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