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_CoupledAggregationFactory.hpp> 55 #include <MueLu_CoupledRBMFactory.hpp> 56 #include <MueLu_DirectSolver.hpp> 57 #include <MueLu_GenericRFactory.hpp> 58 #include <MueLu_Hierarchy.hpp> 59 #include <MueLu_Ifpack2Smoother.hpp> 60 #include <MueLu_PFactory.hpp> 61 #include <MueLu_PgPFactory.hpp> 62 #include <MueLu_RAPFactory.hpp> 63 #include <MueLu_RAPShiftFactory.hpp> 64 #include <MueLu_SaPFactory.hpp> 65 #include <MueLu_ShiftedLaplacian.hpp> 66 #include <MueLu_ShiftedLaplacianOperator.hpp> 67 #include <MueLu_SmootherFactory.hpp> 68 #include <MueLu_SmootherPrototype.hpp> 69 #include <MueLu_TentativePFactory.hpp> 70 #include <MueLu_TransPFactory.hpp> 71 #include <MueLu_UncoupledAggregationFactory.hpp> 72 #include <MueLu_Utilities.hpp> 77 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 coarseGridSize_ = paramList->
get(
"MueLu: coarse size", 1000);
86 numLevels_ = paramList->
get(
"MueLu: levels", 3);
87 int stype = paramList->
get(
"MueLu: smoother", 8);
88 if(stype==1) { Smoother_=
"jacobi"; }
89 else if(stype==2) { Smoother_=
"gauss-seidel"; }
90 else if(stype==3) { Smoother_=
"symmetric gauss-seidel"; }
91 else if(stype==4) { Smoother_=
"chebyshev"; }
92 else if(stype==5) { Smoother_=
"krylov"; }
93 else if(stype==6) { Smoother_=
"ilut"; }
94 else if(stype==7) { Smoother_=
"riluk"; }
95 else if(stype==8) { Smoother_=
"schwarz"; }
96 else if(stype==9) { Smoother_=
"superilu"; }
97 else if(stype==10) { Smoother_=
"superlu"; }
98 else { Smoother_=
"schwarz"; }
99 smoother_sweeps_ = paramList->
get(
"MueLu: sweeps", 5);
100 smoother_damping_ = paramList->
get(
"MueLu: relax val", 1.0);
101 ncycles_ = paramList->
get(
"MueLu: cycles", 1);
102 iters_ = paramList->
get(
"MueLu: iterations", 500);
103 solverType_ = paramList->
get(
"MueLu: solver type", 1);
104 restart_size_ = paramList->
get(
"MueLu: restart size", 100);
105 recycle_size_ = paramList->
get(
"MueLu: recycle size", 25);
106 isSymmetric_ = paramList->
get(
"MueLu: symmetric",
true);
107 ilu_leveloffill_ = paramList->
get(
"MueLu: level-of-fill", 5);
108 ilu_abs_thresh_ = paramList->
get(
"MueLu: abs thresh", 0.0);
109 ilu_rel_thresh_ = paramList->
get(
"MueLu: rel thresh", 1.0);
110 ilu_diagpivotthresh_ = paramList->
get(
"MueLu: piv thresh", 0.1);
111 ilu_drop_tol_ = paramList->
get(
"MueLu: drop tol", 0.01);
112 ilu_fill_tol_ = paramList->
get(
"MueLu: fill tol", 0.01);
113 schwarz_overlap_ = paramList->
get(
"MueLu: overlap", 0);
114 schwarz_usereorder_ = paramList->
get(
"MueLu: use reorder",
true);
115 int combinemode = paramList->
get(
"MueLu: combine mode", 1);
116 if(combinemode==0) { schwarz_combinemode_ = Tpetra::ZERO; }
117 else { schwarz_combinemode_ = Tpetra::ADD; }
118 tol_ = paramList->
get(
"MueLu: tolerance", 0.001);
122 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 if(A_!=Teuchos::null)
128 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT 129 if(LinearProblem_!=Teuchos::null)
130 LinearProblem_ -> setOperator ( TpetraA_ );
137 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
141 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT 142 if(LinearProblem_!=Teuchos::null)
143 LinearProblem_ -> setOperator ( TpetraA_ );
148 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
152 GridTransfersExist_=
false;
156 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
160 =
rcp(
new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraP) );
161 P_=
rcp(
new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp) );
162 GridTransfersExist_=
false;
166 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
173 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
177 =
rcp(
new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraK) );
178 K_=
rcp(
new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp) );
182 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
189 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
193 =
rcp(
new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraM) );
194 M_=
rcp(
new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp) );
198 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
205 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
208 NullSpace_=NullSpace;
212 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
215 levelshifts_=levelshifts;
216 numLevels_=levelshifts_.size();
221 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
235 if(isSymmetric_==
true) {
236 Manager_ -> SetFactory(
"P", Pfact_);
237 Manager_ -> SetFactory(
"R", TransPfact_);
240 Manager_ -> SetFactory(
"P", PgPfact_);
241 Manager_ -> SetFactory(
"R", Rfact_);
244 Manager_ -> SetFactory(
"Ptent", TentPfact_);
246 params.
set(
"lightweight wrap",
true);
247 params.
set(
"aggregation: drop scheme",
"classical");
248 Dropfact_ -> SetParameterList(params);
249 Manager_ -> SetFactory(
"Graph", Dropfact_);
250 if(Aggregation_==
"coupled") {
251 Manager_ -> SetFactory(
"Aggregates", Aggfact_ );
254 Manager_ -> SetFactory(
"Aggregates", UCaggfact_ );
258 if(Smoother_==
"jacobi") {
259 precType_ =
"RELAXATION";
260 precList_.set(
"relaxation: type",
"Jacobi");
261 precList_.set(
"relaxation: sweeps", smoother_sweeps_);
262 precList_.set(
"relaxation: damping factor", smoother_damping_);
264 else if(Smoother_==
"gauss-seidel") {
265 precType_ =
"RELAXATION";
266 precList_.set(
"relaxation: type",
"Gauss-Seidel");
267 precList_.set(
"relaxation: sweeps", smoother_sweeps_);
268 precList_.set(
"relaxation: damping factor", smoother_damping_);
270 else if(Smoother_==
"symmetric gauss-seidel") {
271 precType_ =
"RELAXATION";
272 precList_.set(
"relaxation: type",
"Symmetric Gauss-Seidel");
273 precList_.set(
"relaxation: sweeps", smoother_sweeps_);
274 precList_.set(
"relaxation: damping factor", smoother_damping_);
276 else if(Smoother_==
"chebyshev") {
277 precType_ =
"CHEBYSHEV";
279 else if(Smoother_==
"krylov") {
280 precType_ =
"KRYLOV";
281 precList_.set(
"krylov: iteration type", krylov_type_);
282 precList_.set(
"krylov: number of iterations", krylov_iterations_);
283 precList_.set(
"krylov: residual tolerance",1.0e-8);
284 precList_.set(
"krylov: block size",1);
285 precList_.set(
"krylov: preconditioner type", krylov_preconditioner_);
286 precList_.set(
"relaxation: sweeps",1);
289 else if(Smoother_==
"ilut") {
291 precList_.set(
"fact: ilut level-of-fill", ilu_leveloffill_);
292 precList_.set(
"fact: absolute threshold", ilu_abs_thresh_);
293 precList_.set(
"fact: relative threshold", ilu_rel_thresh_);
294 precList_.set(
"fact: drop tolerance", ilu_drop_tol_);
295 precList_.set(
"fact: relax value", ilu_relax_val_);
297 else if(Smoother_==
"riluk") {
299 precList_.set(
"fact: iluk level-of-fill", ilu_leveloffill_);
300 precList_.set(
"fact: absolute threshold", ilu_abs_thresh_);
301 precList_.set(
"fact: relative threshold", ilu_rel_thresh_);
302 precList_.set(
"fact: drop tolerance", ilu_drop_tol_);
303 precList_.set(
"fact: relax value", ilu_relax_val_);
305 else if(Smoother_==
"schwarz") {
306 precType_ =
"SCHWARZ";
307 precList_.set(
"schwarz: overlap level", schwarz_overlap_);
308 precList_.set(
"schwarz: combine mode", schwarz_combinemode_);
309 precList_.set(
"schwarz: use reordering", schwarz_usereorder_);
311 precList_.set(
"order_method",schwarz_ordermethod_);
312 precList_.sublist(
"schwarz: reordering list").set(
"order_method",schwarz_ordermethod_);
313 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: ilut level-of-fill", ilu_leveloffill_);
314 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: absolute threshold", ilu_abs_thresh_);
315 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: relative threshold", ilu_rel_thresh_);
316 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: drop tolerance", ilu_drop_tol_);
317 precList_.sublist(
"schwarz: subdomain solver parameters").set(
"fact: relax value", ilu_relax_val_);
319 else if(Smoother_==
"superilu") {
320 precType_ =
"superlu";
321 precList_.set(
"RowPerm", ilu_rowperm_);
322 precList_.set(
"ColPerm", ilu_colperm_);
323 precList_.set(
"DiagPivotThresh", ilu_diagpivotthresh_);
324 precList_.set(
"ILU_DropRule",ilu_drop_rule_);
325 precList_.set(
"ILU_DropTol",ilu_drop_tol_);
326 precList_.set(
"ILU_FillFactor",ilu_leveloffill_);
327 precList_.set(
"ILU_Norm",ilu_normtype_);
328 precList_.set(
"ILU_MILU",ilu_milutype_);
329 precList_.set(
"ILU_FillTol",ilu_fill_tol_);
330 precList_.set(
"ILU_Flag",
true);
332 else if(Smoother_==
"superlu") {
333 precType_ =
"superlu";
334 precList_.set(
"ColPerm", ilu_colperm_);
335 precList_.set(
"DiagPivotThresh", ilu_diagpivotthresh_);
337 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT 341 #if defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLU) 342 coarsestSmooProto_ =
rcp(
new DirectSolver(
"Superlu",coarsestSmooList_) );
343 #elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_KLU2) 345 #elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLUDIST) 346 coarsestSmooProto_ =
rcp(
new DirectSolver(
"Superludist",coarsestSmooList_) );
358 if(K_!=Teuchos::null) {
359 Manager_ -> SetFactory(
"Smoother", Teuchos::null);
360 Manager_ -> SetFactory(
"CoarseSolver", Teuchos::null);
362 if(NullSpace_!=Teuchos::null)
363 Hierarchy_ -> GetLevel(0) -> Set(
"Nullspace", NullSpace_);
364 if(isSymmetric_==
true) {
365 Hierarchy_ ->
Keep(
"P", Pfact_.get());
366 Hierarchy_ ->
Keep(
"R", TransPfact_.get());
367 Hierarchy_ -> SetImplicitTranspose(
true);
370 Hierarchy_ ->
Keep(
"P", PgPfact_.get());
371 Hierarchy_ ->
Keep(
"R", Rfact_.get());
373 Hierarchy_ ->
Keep(
"Ptent", TentPfact_.get());
374 Hierarchy_ -> SetMaxCoarseSize( coarseGridSize_ );
375 Hierarchy_ -> Setup(*Manager_, 0, numLevels_);
376 GridTransfersExist_=
true;
380 Manager_ -> SetFactory(
"Smoother", smooFact_);
381 Manager_ -> SetFactory(
"CoarseSolver", coarsestSmooFact_);
383 if(NullSpace_!=Teuchos::null)
384 Hierarchy_ -> GetLevel(0) -> Set(
"Nullspace", NullSpace_);
385 if(isSymmetric_==
true)
386 Hierarchy_ -> SetImplicitTranspose(
true);
387 Hierarchy_ -> SetMaxCoarseSize( coarseGridSize_ );
388 Hierarchy_ -> Setup(*Manager_, 0, numLevels_);
389 GridTransfersExist_=
true;
394 BelosList_ ->
set(
"Maximum Iterations",iters_ );
395 BelosList_ ->
set(
"Convergence Tolerance",tol_ );
397 BelosList_ ->
set(
"Output Frequency",1);
398 BelosList_ ->
set(
"Output Style",Belos::Brief);
399 BelosList_ ->
set(
"Num Blocks",restart_size_);
400 BelosList_ ->
set(
"Num Recycled Blocks",recycle_size_);
407 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
410 int numLevels = Hierarchy_ -> GetNumLevels();
411 Manager_ -> SetFactory(
"Smoother", smooFact_);
412 Manager_ -> SetFactory(
"CoarseSolver", coarsestSmooFact_);
413 Hierarchy_ -> GetLevel(0) -> Set(
"A", P_);
414 Hierarchy_ -> Setup(*Manager_, 0, numLevels);
420 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
423 int numLevels = Hierarchy_ -> GetNumLevels();
424 Acshift_->SetShifts(levelshifts_);
425 Manager_ -> SetFactory(
"Smoother", smooFact_);
426 Manager_ -> SetFactory(
"CoarseSolver", coarsestSmooFact_);
427 Manager_ -> SetFactory(
"A", Acshift_);
428 Manager_ -> SetFactory(
"K", Acshift_);
429 Manager_ -> SetFactory(
"M", Acshift_);
430 Hierarchy_ -> GetLevel(0) -> Set(
"A", P_);
431 Hierarchy_ -> GetLevel(0) -> Set(
"K", K_);
432 Hierarchy_ -> GetLevel(0) -> Set(
"M", M_);
433 Hierarchy_ -> Setup(*Manager_, 0, numLevels);
439 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
443 if( GridTransfersExist_ ==
false ) {
445 if(NullSpace_!=Teuchos::null)
446 Hierarchy_ -> GetLevel(0) -> Set(
"Nullspace", NullSpace_);
447 if(isSymmetric_==
true)
448 Hierarchy_ -> SetImplicitTranspose(
true);
449 Hierarchy_ -> SetMaxCoarseSize( coarseGridSize_ );
450 Hierarchy_ -> Setup(*Manager_, 0, numLevels_);
451 GridTransfersExist_=
true;
457 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
460 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT 463 (Hierarchy_, A_, ncycles_, subiters_, option_, tol_) );
465 if(LinearProblem_==Teuchos::null)
466 LinearProblem_ =
rcp(
new LinearProblem );
467 LinearProblem_ -> setOperator ( TpetraA_ );
468 LinearProblem_ -> setRightPrec( MueLuOp_ );
469 if(SolverManager_==Teuchos::null) {
470 std::string solverName;
471 SolverFactory_=
rcp(
new SolverFactory() );
472 if(solverType_==1) { solverName=
"Block GMRES"; }
473 else if(solverType_==2) { solverName=
"Recycling GMRES"; }
474 else { solverName=
"Flexible GMRES"; }
475 SolverManager_ = SolverFactory_->create( solverName, BelosList_ );
476 SolverManager_ -> setProblem( LinearProblem_ );
483 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
486 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT 487 LinearProblem_ -> setOperator ( TpetraA_ );
494 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
497 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT 499 LinearProblem_ -> setProblem(X, B);
501 SolverManager_ -> solve();
509 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
514 Hierarchy_ -> Iterate(*B, *X, 1,
true, 0);
518 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
520 RCP<Tpetra::MultiVector<SC,LO,GO,NO> >& X)
523 =
Teuchos::rcp(
new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(X) );
525 =
Teuchos::rcp(
new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(B) );
527 Hierarchy_ -> Iterate(*XpetraB, *XpetraX, 1,
true, 0);
531 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
534 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT 535 int numiters = SolverManager_ -> getNumIters();
544 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
547 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT 548 double residual = SolverManager_ -> achievedTol();
558 #define MUELU_SHIFTEDLAPLACIAN_SHORT 560 #endif //if defined(HAVE_MUELU_IFPACK2) and defined(HAVE_MUELU_TPETRA) 561 #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 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)
Namespace for MueLu classes and methods.
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_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.Input/output of SaPFactory
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)