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_CoalesceDropFactory.hpp>
54 #include <MueLu_CoarseMapFactory.hpp>
55 #include <MueLu_CoupledAggregationFactory.hpp>
56 #include <MueLu_CoupledRBMFactory.hpp>
57 #include <MueLu_DirectSolver.hpp>
58 #include <MueLu_GenericRFactory.hpp>
59 #include <MueLu_Hierarchy.hpp>
60 #include <MueLu_Ifpack2Smoother.hpp>
61 #include <MueLu_PFactory.hpp>
62 #include <MueLu_PgPFactory.hpp>
63 #include <MueLu_RAPFactory.hpp>
64 #include <MueLu_RAPShiftFactory.hpp>
65 #include <MueLu_SaPFactory.hpp>
66 #include <MueLu_ShiftedLaplacian.hpp>
67 #include <MueLu_ShiftedLaplacianOperator.hpp>
68 #include <MueLu_SmootherFactory.hpp>
69 #include <MueLu_SmootherPrototype.hpp>
70 #include <MueLu_TentativePFactory.hpp>
71 #include <MueLu_TransPFactory.hpp>
72 #include <MueLu_UncoupledAggregationFactory.hpp>
73 #include <MueLu_Utilities.hpp>
78 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
82 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
123 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
127 if(A_!=Teuchos::null)
129 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
130 if(LinearProblem_!=Teuchos::null)
131 LinearProblem_ -> setOperator ( TpetraA_ );
138 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
142 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
143 if(LinearProblem_!=Teuchos::null)
144 LinearProblem_ -> setOperator ( TpetraA_ );
149 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
153 GridTransfersExist_=
false;
157 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
163 GridTransfersExist_=
false;
167 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
174 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
183 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
190 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
199 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
206 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
209 NullSpace_=NullSpace;
213 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
216 levelshifts_=levelshifts;
217 numLevels_=levelshifts_.size();
222 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
237 if(isSymmetric_==
true) {
238 Manager_ -> SetFactory(
"P", Pfact_);
239 Manager_ -> SetFactory(
"R", TransPfact_);
242 Manager_ -> SetFactory(
"P", PgPfact_);
243 Manager_ -> SetFactory(
"R", Rfact_);
246 Manager_ -> SetFactory(
"Ptent", TentPfact_);
248 params.
set(
"lightweight wrap",
true);
249 params.
set(
"aggregation: drop scheme",
"classical");
250 Dropfact_ -> SetParameterList(params);
251 Manager_ -> SetFactory(
"Graph", Dropfact_);
252 if(Aggregation_==
"coupled") {
253 Manager_ -> SetFactory(
"Aggregates", Aggfact_ );
256 Manager_ -> SetFactory(
"Aggregates", UCaggfact_ );
258 Manager_ -> SetFactory(
"CoarseMap", CoarseMapfact_);
261 if(Smoother_==
"jacobi") {
262 precType_ =
"RELAXATION";
263 precList_.set(
"relaxation: type",
"Jacobi");
264 precList_.set(
"relaxation: sweeps", smoother_sweeps_);
265 precList_.set(
"relaxation: damping factor", smoother_damping_);
267 else if(Smoother_==
"gauss-seidel") {
268 precType_ =
"RELAXATION";
269 precList_.set(
"relaxation: type",
"Gauss-Seidel");
270 precList_.set(
"relaxation: sweeps", smoother_sweeps_);
271 precList_.set(
"relaxation: damping factor", smoother_damping_);
273 else if(Smoother_==
"symmetric gauss-seidel") {
274 precType_ =
"RELAXATION";
275 precList_.set(
"relaxation: type",
"Symmetric Gauss-Seidel");
276 precList_.set(
"relaxation: sweeps", smoother_sweeps_);
277 precList_.set(
"relaxation: damping factor", smoother_damping_);
279 else if(Smoother_==
"chebyshev") {
280 precType_ =
"CHEBYSHEV";
282 else if(Smoother_==
"krylov") {
283 precType_ =
"KRYLOV";
284 precList_.set(
"krylov: iteration type", krylov_type_);
285 precList_.set(
"krylov: number of iterations", krylov_iterations_);
286 precList_.set(
"krylov: residual tolerance",1.0e-8);
287 precList_.set(
"krylov: block size",1);
288 precList_.set(
"krylov: preconditioner type", krylov_preconditioner_);
289 precList_.set(
"relaxation: sweeps",1);
292 else if(Smoother_==
"ilut") {
294 precList_.set(
"fact: ilut level-of-fill", ilu_leveloffill_);
295 precList_.set(
"fact: absolute threshold", ilu_abs_thresh_);
296 precList_.set(
"fact: relative threshold", ilu_rel_thresh_);
297 precList_.set(
"fact: drop tolerance", ilu_drop_tol_);
298 precList_.set(
"fact: relax value", ilu_relax_val_);
300 else if(Smoother_==
"riluk") {
302 precList_.set(
"fact: iluk level-of-fill", ilu_leveloffill_);
303 precList_.set(
"fact: absolute threshold", ilu_abs_thresh_);
304 precList_.set(
"fact: relative threshold", ilu_rel_thresh_);
305 precList_.set(
"fact: drop tolerance", ilu_drop_tol_);
306 precList_.set(
"fact: relax value", ilu_relax_val_);
308 else if(Smoother_==
"schwarz") {
309 precType_ =
"SCHWARZ";
310 precList_.set(
"schwarz: overlap level", schwarz_overlap_);
311 precList_.set(
"schwarz: combine mode", schwarz_combinemode_);
312 precList_.set(
"schwarz: use reordering", schwarz_usereorder_);
314 precList_.set(
"order_method",schwarz_ordermethod_);
315 precList_.sublist(
"schwarz: reordering list").set(
"order_method",schwarz_ordermethod_);
316 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: ilut level-of-fill", ilu_leveloffill_);
317 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: absolute threshold", ilu_abs_thresh_);
318 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: relative threshold", ilu_rel_thresh_);
319 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: drop tolerance", ilu_drop_tol_);
320 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: relax value", ilu_relax_val_);
322 else if(Smoother_==
"superilu") {
323 precType_ =
"superlu";
324 precList_.set(
"RowPerm", ilu_rowperm_);
325 precList_.set(
"ColPerm", ilu_colperm_);
326 precList_.set(
"DiagPivotThresh", ilu_diagpivotthresh_);
327 precList_.set(
"ILU_DropRule",ilu_drop_rule_);
328 precList_.set(
"ILU_DropTol",ilu_drop_tol_);
329 precList_.set(
"ILU_FillFactor",ilu_leveloffill_);
330 precList_.set(
"ILU_Norm",ilu_normtype_);
331 precList_.set(
"ILU_MILU",ilu_milutype_);
332 precList_.set(
"ILU_FillTol",ilu_fill_tol_);
333 precList_.set(
"ILU_Flag",
true);
335 else if(Smoother_==
"superlu") {
336 precType_ =
"superlu";
337 precList_.set(
"ColPerm", ilu_colperm_);
338 precList_.set(
"DiagPivotThresh", ilu_diagpivotthresh_);
340 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
344 #if defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLU)
345 coarsestSmooProto_ =
rcp(
new DirectSolver(
"Superlu",coarsestSmooList_) );
346 #elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_KLU2)
348 #elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLUDIST)
349 coarsestSmooProto_ =
rcp(
new DirectSolver(
"Superludist",coarsestSmooList_) );
361 if(K_!=Teuchos::null) {
362 Manager_ -> SetFactory(
"Smoother", Teuchos::null);
363 Manager_ -> SetFactory(
"CoarseSolver", Teuchos::null);
365 if(NullSpace_!=Teuchos::null)
366 Hierarchy_ -> GetLevel(0) -> Set(
"Nullspace", NullSpace_);
367 if(isSymmetric_==
true) {
368 Hierarchy_ ->
Keep(
"P", Pfact_.get());
369 Hierarchy_ ->
Keep(
"R", TransPfact_.get());
370 Hierarchy_ -> SetImplicitTranspose(
true);
373 Hierarchy_ ->
Keep(
"P", PgPfact_.get());
374 Hierarchy_ ->
Keep(
"R", Rfact_.get());
376 Hierarchy_ ->
Keep(
"Ptent", TentPfact_.get());
377 Hierarchy_ -> SetMaxCoarseSize( coarseGridSize_ );
378 Hierarchy_ -> Setup(*Manager_, 0, numLevels_);
379 GridTransfersExist_=
true;
383 Manager_ -> SetFactory(
"Smoother", smooFact_);
384 Manager_ -> SetFactory(
"CoarseSolver", coarsestSmooFact_);
386 if(NullSpace_!=Teuchos::null)
387 Hierarchy_ -> GetLevel(0) -> Set(
"Nullspace", NullSpace_);
388 if(isSymmetric_==
true)
389 Hierarchy_ -> SetImplicitTranspose(
true);
390 Hierarchy_ -> SetMaxCoarseSize( coarseGridSize_ );
391 Hierarchy_ -> Setup(*Manager_, 0, numLevels_);
392 GridTransfersExist_=
true;
397 BelosList_ -> set(
"Maximum Iterations",iters_ );
398 BelosList_ -> set(
"Convergence Tolerance",tol_ );
400 BelosList_ -> set(
"Output Frequency",1);
401 BelosList_ -> set(
"Output Style",Belos::Brief);
402 BelosList_ -> set(
"Num Blocks",restart_size_);
403 BelosList_ -> set(
"Num Recycled Blocks",recycle_size_);
410 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
413 int numLevels = Hierarchy_ -> GetNumLevels();
414 Manager_ -> SetFactory(
"Smoother", smooFact_);
415 Manager_ -> SetFactory(
"CoarseSolver", coarsestSmooFact_);
416 Hierarchy_ -> GetLevel(0) -> Set(
"A", P_);
417 Hierarchy_ -> Setup(*Manager_, 0, numLevels);
423 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
426 int numLevels = Hierarchy_ -> GetNumLevels();
427 Acshift_->SetShifts(levelshifts_);
428 Manager_ -> SetFactory(
"Smoother", smooFact_);
429 Manager_ -> SetFactory(
"CoarseSolver", coarsestSmooFact_);
430 Manager_ -> SetFactory(
"A", Acshift_);
431 Manager_ -> SetFactory(
"K", Acshift_);
432 Manager_ -> SetFactory(
"M", Acshift_);
433 Hierarchy_ -> GetLevel(0) -> Set(
"A", P_);
434 Hierarchy_ -> GetLevel(0) -> Set(
"K", K_);
435 Hierarchy_ -> GetLevel(0) -> Set(
"M", M_);
436 Hierarchy_ -> Setup(*Manager_, 0, numLevels);
442 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
446 if( GridTransfersExist_ ==
false ) {
448 if(NullSpace_!=Teuchos::null)
449 Hierarchy_ -> GetLevel(0) -> Set(
"Nullspace", NullSpace_);
450 if(isSymmetric_==
true)
451 Hierarchy_ -> SetImplicitTranspose(
true);
452 Hierarchy_ -> SetMaxCoarseSize( coarseGridSize_ );
453 Hierarchy_ -> Setup(*Manager_, 0, numLevels_);
454 GridTransfersExist_=
true;
460 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
463 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
466 (Hierarchy_, A_, ncycles_, subiters_, option_, tol_) );
468 if(LinearProblem_==Teuchos::null)
469 LinearProblem_ =
rcp(
new LinearProblem );
470 LinearProblem_ -> setOperator ( TpetraA_ );
471 LinearProblem_ -> setRightPrec( MueLuOp_ );
472 if(SolverManager_==Teuchos::null) {
473 std::string solverName;
474 SolverFactory_=
rcp(
new SolverFactory() );
475 if(solverType_==1) { solverName=
"Block GMRES"; }
476 else if(solverType_==2) { solverName=
"Recycling GMRES"; }
477 else { solverName=
"Flexible GMRES"; }
478 SolverManager_ = SolverFactory_->create( solverName, BelosList_ );
479 SolverManager_ -> setProblem( LinearProblem_ );
486 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
489 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
490 LinearProblem_ -> setOperator ( TpetraA_ );
497 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
500 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
502 LinearProblem_ -> setProblem(X, B);
504 SolverManager_ -> solve();
512 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
517 Hierarchy_ -> Iterate(*B, *X, 1,
true, 0);
521 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
523 RCP<Tpetra::MultiVector<SC,LO,GO,NO> >& X)
530 Hierarchy_ -> Iterate(*XpetraB, *XpetraX, 1,
true, 0);
534 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
537 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
538 int numiters = SolverManager_ -> getNumIters();
547 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
552 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
553 MT residual = SolverManager_ -> achievedTol();
563 #define MUELU_SHIFTEDLAPLACIAN_SHORT
565 #endif //if defined(HAVE_MUELU_IFPACK2) and defined(HAVE_MUELU_TPETRA)
566 #endif // MUELU_SHIFTEDLAPLACIAN_DEF_HPP
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
This class specifies the default factory that should generate some data on a Level if the data does n...
Factory for generating coarse level map. Used by TentativePFactory.
Factory for building coarse grid matrices, when the matrix is of the form K+a*M. Useful when you want...
void multigrid_apply(const RCP< MultiVector > B, RCP< MultiVector > &X)
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setmass(RCP< Matrix > &M)
void resetLinearProblem()
Factory for building tentative prolongator.
void setParameters(Teuchos::RCP< Teuchos::ParameterList > paramList)
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
Factory for coarsening a graph with uncoupled aggregation.
Factory for building restriction operators using a prolongator factory.
Teuchos::ScalarTraits< Scalar >::magnitudeType GetResidual()
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setLevelShifts(std::vector< Scalar > levelshifts)
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
virtual ~ShiftedLaplacian()
Always keep data, even accross run. This flag is set by Level::Keep(). This flag is propagated to coa...
void setcoords(RCP< MultiVector > &Coords)
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator, with an optional two-level correction...
int solve(const RCP< TMV > B, RCP< TMV > &X)
Factory for creating a graph base on a given matrix.
void setProblemMatrix(RCP< Matrix > &A)
void setstiff(RCP< Matrix > &K)
Factory for building restriction operators.
Exception throws to report errors in the internal logical of the program.
Print all warning messages.
Factory for building coarse matrices.
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
Class that encapsulates Ifpack2 smoothers.
Factory for building Smoothed Aggregation prolongators.
Factory for building uncoupled aggregates.
void setNullSpace(RCP< MultiVector > NullSpace)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
void setPreconditioningMatrix(RCP< Matrix > &P)