46#ifndef MUELU_SHIFTEDLAPLACIAN_DEF_HPP
47#define MUELU_SHIFTEDLAPLACIAN_DEF_HPP
51#if defined(HAVE_MUELU_IFPACK2) and defined(HAVE_MUELU_TPETRA)
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>
79template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
124template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
128 if(A_!=Teuchos::null)
130#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
131 if(LinearProblem_!=Teuchos::null)
132 LinearProblem_ -> setOperator ( TpetraA_ );
134 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"ShiftedLaplacian only available with Tpetra and GO=int enabled.");
139template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
143#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
144 if(LinearProblem_!=Teuchos::null)
145 LinearProblem_ -> setOperator ( TpetraA_ );
150template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
154 GridTransfersExist_=
false;
158template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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;
168template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
175template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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) );
184template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
191template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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) );
200template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
207template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
210 NullSpace_=NullSpace;
214template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
217 levelshifts_=levelshifts;
218 numLevels_=levelshifts_.size();
223template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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_ );
249 Manager_ -> SetFactory(
"Aggregates", UCaggfact_ );
251 Manager_ -> SetFactory(
"CoarseMap", CoarseMapfact_);
252 Manager_ -> SetFactory(
"Ptent", TentPfact_);
253 if(isSymmetric_==
true) {
254 Manager_ -> SetFactory(
"P", Pfact_);
255 Manager_ -> SetFactory(
"R", TransPfact_);
258 Manager_ -> SetFactory(
"P", PgPfact_);
259 Manager_ -> SetFactory(
"R", Rfact_);
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_);
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_);
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_);
282 else if(Smoother_==
"chebyshev") {
283 precType_ =
"CHEBYSHEV";
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);
295 else if(Smoother_==
"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_);
303 else if(Smoother_==
"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_);
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_);
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_);
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);
338 else if(Smoother_==
"superlu") {
339 precType_ =
"superlu";
340 precList_.set(
"ColPerm", ilu_colperm_);
341 precList_.set(
"DiagPivotThresh", ilu_diagpivotthresh_);
343#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
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_) );
356 coarsestSmooFact_ = rcp(
new SmootherFactory(coarsestSmooProto_, Teuchos::null) );
364 if(K_!=Teuchos::null) {
365 Manager_ -> SetFactory(
"Smoother", Teuchos::null);
366 Manager_ -> SetFactory(
"CoarseSolver", Teuchos::null);
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);
376 Hierarchy_ ->
Keep(
"P", PgPfact_.get());
377 Hierarchy_ ->
Keep(
"R", Rfact_.get());
379 Hierarchy_ ->
Keep(
"Ptent", TentPfact_.get());
380 Hierarchy_ -> SetMaxCoarseSize( coarseGridSize_ );
381 Hierarchy_ -> Setup(*Manager_, 0, numLevels_);
382 GridTransfersExist_=
true;
386 Manager_ -> SetFactory(
"Smoother", smooFact_);
387 Manager_ -> SetFactory(
"CoarseSolver", coarsestSmooFact_);
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;
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_);
408 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"ShiftedLaplacian only available with Tpetra and GO=int enabled.");
413template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
426template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
445template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
449 if( GridTransfersExist_ ==
false ) {
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;
463template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
466#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
469 (Hierarchy_, A_, ncycles_, subiters_, option_, tol_) );
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_ );
485 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"ShiftedLaplacian only available with Tpetra and GO=int enabled.");
489template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
492#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
493 LinearProblem_ -> setOperator ( TpetraA_ );
495 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"ShiftedLaplacian only available with Tpetra and GO=int enabled.");
500template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
503#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
505 LinearProblem_ -> setProblem(X, B);
507 SolverManager_ -> solve();
509 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"ShiftedLaplacian only available with Tpetra and GO=int enabled.");
515template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
520 Hierarchy_ -> Iterate(*B, *X, 1,
true, 0);
524template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
526 RCP<Tpetra::MultiVector<SC,LO,GO,NO> >& X)
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) );
533 Hierarchy_ -> Iterate(*XpetraB, *XpetraX, 1,
true, 0);
537template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
540#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
541 int numiters = SolverManager_ -> getNumIters();
544 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"ShiftedLaplacian only available with Tpetra and GO=int enabled.");
550template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
551typename Teuchos::ScalarTraits<Scalar>::magnitudeType
554 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
555#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
556 MT residual = SolverManager_ -> achievedTol();
559 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"ShiftedLaplacian only available with Tpetra and GO=int enabled.");
566#define MUELU_SHIFTEDLAPLACIAN_SHORT
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 setmass(RCP< Matrix > &M)
void resetLinearProblem()
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)
virtual ~ShiftedLaplacian()
void setstiff(RCP< Matrix > &K)
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.
Factory for building uncoupled aggregates.
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...