46 #ifndef MUELU_REFMAXWELL_DEF_HPP 47 #define MUELU_REFMAXWELL_DEF_HPP 51 #include "Xpetra_MatrixMatrix.hpp" 53 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2) 56 #include "MueLu_Utilities.hpp" 58 #include "MueLu_MLParameterListInterpreter.hpp" 59 #include "MueLu_ParameterListInterpreter.hpp" 64 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67 return Xpetra::toTpetraNonZero(SM_Matrix_->getDomainMap());
71 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
74 return Xpetra::toTpetraNonZero(SM_Matrix_->getRangeMap());
78 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 disable_addon_ = list.
get(
"refmaxwell: disable add-on",
true);
82 mode_ = list.
get(
"refmaxwell: mode",
"additive");
85 precList11_ = list.
sublist(
"refmaxwell: 11list");
88 precList22_ = list.
sublist(
"refmaxwell: 22list");
90 std::string ref(
"smoother:");
91 std::string replace(
"coarse:");
93 const std::string & pname = list.
name(i);
94 if(pname.find(ref)!=std::string::npos) {
95 smootherList_.setEntry(pname,list.
entry(i));
96 std::string coarsename(pname);
97 coarsename.replace((
size_t)0,(
size_t)ref.length(),replace);
101 smootherList_.set(
"coarse: params",list.
sublist(
"smoother: params"));
103 smootherList_.set(
"max levels",1);
107 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
115 findDirichletRows(SM_Matrix_,BCrows_);
116 findDirichletCols(D0_Matrix_,BCrows_,BCcols_);
117 D0_Matrix_->resumeFill();
118 Apply_BCsToMatrixRows(D0_Matrix_,BCrows_);
119 Apply_BCsToMatrixCols(D0_Matrix_,BCcols_);
120 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
125 TMT_Matrix_=MatrixFactory::Build(D0_Matrix_->getDomainMap(),0);
126 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*SM_Matrix_,
false,*D0_Matrix_,
false,*C1,
true,
true);
127 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
true,*C1,
false,*TMT_Matrix_,
true,
true);
128 TMT_Matrix_->resumeFill();
129 Remove_Zeroed_Rows(TMT_Matrix_,1.0e-16);
130 TMT_Matrix_->SetFixedBlockSize(1);
134 if(Nullspace_ != Teuchos::null) {
137 else if(Nullspace_ == Teuchos::null && Coords_ != Teuchos::null) {
138 Nullspace_ = MultiVectorFactory::Build(SM_Matrix_->getRowMap(),Coords_->getNumVectors());
139 D0_Matrix_->apply(*Coords_,*Nullspace_);
142 std::cerr <<
"MueLu::RefMaxwell::compute(): either the nullspace or the nodal coordinates must be provided." << std::endl;
146 if(P11_==Teuchos::null) {
155 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M1_Matrix_,
false,*D0_Matrix_,
false,*C,
true,
true);
156 A22_=MatrixFactory::Build(D0_Matrix_->getDomainMap(),0);
157 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
true,*C,
false,*A22_,
true,
true);
159 Remove_Zeroed_Rows(A22_,1.0e-16);
160 A22_->SetFixedBlockSize(1);
165 std::string syntaxStr =
"parameterlist: syntax";
166 if (parameterList_.isParameter(syntaxStr) && parameterList_.get<std::string>(syntaxStr) ==
"ml") {
167 parameterList_.remove(syntaxStr);
178 Hierarchy11_->setlib(Xpetra::UseTpetra);
179 Hierarchy11_->GetLevel(0)->Set(
"A", A11_);
183 Hierarchy22_->setlib(Xpetra::UseTpetra);
184 Hierarchy22_->GetLevel(0)->Set(
"A", A22_);
189 HierarchySmoother_->setlib(Xpetra::UseTpetra);
190 HierarchySmoother_->GetLevel(0)->Set(
"A", SM_Matrix_);
195 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
215 params.
set(
"sa: damping factor",0.0);
216 Pfact -> SetParameterList(params);
217 auxManager -> SetFactory(
"P", Pfact);
218 auxManager -> SetFactory(
"Ptent", TentPfact);
219 auxManager -> SetFactory(
"Aggregates", Aggfact);
220 auxManager -> SetFactory(
"Smoother", Teuchos::null);
221 auxManager -> SetFactory(
"CoarseSolver", Teuchos::null);
222 auxHierarchy ->
Keep(
"P", Pfact.
get());
223 auxHierarchy -> SetMaxCoarseSize(1);
224 auxHierarchy -> Setup(*auxManager, 0, 2);
232 D0_Matrix_Abs -> resumeFill();
233 D0_Matrix_Abs -> setAllToScalar((Scalar)0.5);
234 Apply_BCsToMatrixRows(D0_Matrix_Abs,BCrows_);
235 Apply_BCsToMatrixCols(D0_Matrix_Abs,BCcols_);
236 D0_Matrix_Abs -> fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
238 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_Abs,
false,*P,
false,*Ptent,
true,
true);
241 size_t dim = Nullspace_->getNumVectors();
242 size_t numLocalRows = SM_Matrix_->getNodeNumRows();
244 = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(Ptent->getColMap(),dim);
247 std::vector< Teuchos::ArrayRCP<const Scalar> > nullspace(dim);
248 for(
size_t i=0; i<dim; i++) {
250 nullspace[i]=datavec;
254 std::vector<size_t> rowPtrs;
255 std::vector<LocalOrdinal> blockCols;
256 std::vector<Scalar> blockVals;
257 for(
size_t i=0; i<numLocalRows; i++) {
258 rowPtrs.push_back(nnz);
261 Ptent->getLocalRowView(i,localCols,localVals);
262 size_t numCols = localCols.
size();
263 for(
size_t j=0; j<numCols; j++) {
264 for(
size_t k=0; k<dim; k++) {
265 blockCols.push_back(localCols[j]*dim+k);
266 blockVals.push_back(localVals[j]*nullspace[k][i]);
271 rowPtrs.push_back(nnz);
278 TP11->allocateAllValues(nnz, rcpRowPtr, rcpColumns, rcpValues);
284 for (
size_t ii = 0; ii < rowPtrs.size(); ii++) rows[ii] = rowPtrs[ii];
285 for (
size_t ii = 0; ii < blockCols.size(); ii++) columns[ii] = blockCols[ii];
286 for (
size_t ii = 0; ii < blockVals.size(); ii++) values[ii] = blockVals[ii];
287 TP11->setAllValues(rcpRowPtr, rcpColumns, rcpValues);
289 = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(Ptent->getDomainMap(),dim);
290 TP11->expertStaticFillComplete(blockCoarseMap,SM_Matrix_->getDomainMap());
294 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
302 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*SM_Matrix_,
false,*P11_,
false,*C,
true,
true);
305 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*P11_,
true,*C,
false,*Matrix1,
true,
true);
307 if(disable_addon_==
true) {
314 "MueLu::RefMaxwell::formCoarseMatrix(): Inverse of " 315 "lumped mass matrix required for add-on (i.e. M0inv_Matrix is null)");
322 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M1_Matrix_,
false,*P11_,
false,*Zaux,
true,
true);
324 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
true,*Zaux,
false,*Z,
true,
true);
326 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M0inv_Matrix_,
false,*Z,
false,*C2,
true,
true);
329 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*Z,
true,*C2,
false,*Matrix2,
true,
true);
332 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*Matrix1,
false,(Scalar)1.0,*Matrix2,
false,(Scalar)1.0,A11_,*out);
333 A11_->fillComplete();
337 size_t dim = Nullspace_->getNumVectors();
338 A11_->SetFixedBlockSize(dim);
342 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
351 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
356 RCP<XMV> P11res = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
357 RCP<XMV> P11x = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
358 RCP<XMV> D0res = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
359 RCP<XMV> D0x = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
364 Hierarchy11_->Iterate(*P11res, *P11x, 1,
true);
365 Hierarchy22_->Iterate(*D0res, *D0x, 1,
true);
370 X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
374 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
377 RCP<XMV> P11res = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
378 RCP<XMV> P11x = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
379 RCP<XMV> D0res = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
380 RCP<XMV> D0x = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
385 Hierarchy11_->Iterate(*P11res, *P11x, 1,
true);
387 X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
392 Hierarchy22_->Iterate(*D0res, *D0x, 1,
true);
394 X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
399 Hierarchy11_->Iterate(*P11res, *P11x, 1,
true);
401 X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
405 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
408 RCP<XMV> P11res = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
409 RCP<XMV> P11x = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
410 RCP<XMV> D0res = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
411 RCP<XMV> D0x = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
416 Hierarchy22_->Iterate(*D0res, *D0x, 1,
true);
418 X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
423 Hierarchy11_->Iterate(*P11res, *P11x, 1,
true);
425 X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
430 Hierarchy22_->Iterate(*D0res, *D0x, 1,
true);
432 X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
436 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
438 Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
443 TMV& temp_x =
const_cast<TMV &
>(X);
444 const XTMV tX(rcpFromRef(temp_x));
445 XTMV tY(rcpFromRef(Y));
449 HierarchySmoother_->Iterate(tX,tY,1);
452 if(mode_==
"additive")
453 applyInverseAdditive(tX,tY);
454 else if(mode_==
"121")
455 applyInverse121(tX,tY);
456 else if(mode_==
"212")
457 applyInverse212(tX,tY);
459 applyInverseAdditive(tX,tY);
462 HierarchySmoother_->Iterate(tX,tY,1);
464 }
catch (std::exception& e) {
467 std::cerr <<
"Caught an exception in MueLu::RefMaxwell::ApplyInverse():" << std::endl
468 << e.what() << std::endl;
473 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
478 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
491 Hierarchy11_ = Teuchos::null;
492 Hierarchy22_ = Teuchos::null;
493 HierarchySmoother_ = Teuchos::null;
494 parameterList_ = List;
495 disable_addon_ =
false;
505 if(M0inv_Matrix != Teuchos::null) {
514 if(Coords != Teuchos::null)
515 Coords_ = Xpetra::toXpetra(Coords);
516 if(Nullspace != Teuchos::null)
517 Nullspace_ = Xpetra::toXpetra(Nullspace);
521 #endif //ifdef HAVE_MUELU_TPETRA 523 #define MUELU_REFMAXWELL_SHORT 524 #endif //ifdef MUELU_REFMAXWELL_DEF_HPP This class specifies the default factory that should generate some data on a Level if the data does n...
basic_FancyOStream & setShowProcRank(const bool showProcRank)
void applyInverse121(const XTMV &RHS, XTMV &X) const
apply 1-2-1 algorithm for 2x2 solve
ConstIterator begin() const
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)
const ParameterEntry & entry(ConstIterator i) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool isSublist(const std::string &name) const
Namespace for MueLu classes and methods.
void setParameters(Teuchos::ParameterList &list)
Set parameters.
void buildProlongator()
Setup the prolongator for the (1,1)-block.
Factory for building tentative prolongator.
Xpetra::TpetraMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > XTMV
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.
void applyInverse212(const XTMV &RHS, XTMV &X) const
apply 2-1-2 algorithm for 2x2 solve
ConstIterator end() const
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra::CrsMatrixWrap< Scalar, LocalOrdinal, GlobalOrdinal, Node > XCrsWrap
void compute()
Setup the preconditioner.
params_t::ConstIterator ConstIterator
Always keep data, even accross run. This flag is set by Level::Keep(). This flag is propagated to coa...
void applyInverseAdditive(const XTMV &RHS, XTMV &X) const
apply additive algorithm for 2x2 solve
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
const std::string & name() const
Class that accepts ML-style parameters and builds a MueLu preconditioner. This interpreter uses the s...
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
virtual RCP< Hierarchy > CreateHierarchy() const
Create an empty Hierarchy object.
Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > TMV
void initialize(const Teuchos::RCP< TCRS > &D0_Matrix, const Teuchos::RCP< TCRS > &M0inv_Matrix, const Teuchos::RCP< TCRS > &M1_Matrix, const Teuchos::RCP< TMV > &Nullspace, const Teuchos::RCP< TMV > &Coords, Teuchos::ParameterList &List)
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Xpetra::TpetraCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > XTCRS
void apply(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::one()) const
void formCoarseMatrix()
Compute P11^{T}*A*P11 efficiently.
#define TEUCHOS_ASSERT(assertion_test)
void resetMatrix(Teuchos::RCP< TCRS > SM_Matrix_new)
Reset system matrix.
Factory for building Smoothed Aggregation prolongators.Input/output of SaPFactory
Factory for building uncoupled aggregates.
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.