46 #ifndef MUELU_SAPFACTORY_DEF_HPP 47 #define MUELU_SAPFACTORY_DEF_HPP 49 #include <Xpetra_Matrix.hpp> 50 #include <Xpetra_IteratorOps.hpp> 51 #include <Xpetra_IO.hpp> 60 #include "MueLu_PerfUtils.hpp" 62 #include "MueLu_TentativePFactory.hpp" 63 #include "MueLu_Utilities.hpp" 67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
69 RCP<ParameterList> validParamList = rcp(
new ParameterList());
71 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 81 #undef SET_VALID_ENTRY 83 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
84 validParamList->set< RCP<const FactoryBase> >(
"P", Teuchos::null,
"Tentative prolongator factory");
87 ParameterList norecurse;
88 norecurse.disableRecursiveValidation();
89 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
92 return validParamList;
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
97 Input(fineLevel,
"A");
101 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
102 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.
GetFactoryManager()->GetFactory(
"Ptent"); }
103 coarseLevel.
DeclareInput(
"P", initialPFact.get(),
this);
106 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
108 return BuildP(fineLevel, coarseLevel);
111 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 const std::string prefix =
"MueLu::SaPFactory(" + levelIDs +
"): ";
119 typedef typename Teuchos::ScalarTraits<SC>::coordinateType Coordinate;
124 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
125 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.
GetFactoryManager()->GetFactory(
"Ptent"); }
126 const ParameterList& pL = GetParameterList();
129 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
130 RCP<Matrix> Ptent = coarseLevel.
Get< RCP<Matrix> >(
"P", initialPFact.get());
132 if (restrictionMode_) {
141 RCP<ParameterList> APparams = rcp(
new ParameterList);;
142 if(pL.isSublist(
"matrixmatrix: kernel params"))
143 APparams->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
145 if (coarseLevel.
IsAvailable(
"AP reuse data",
this)) {
146 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
148 APparams = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data",
this);
150 if (APparams->isParameter(
"graph"))
151 finalP = APparams->get< RCP<Matrix> >(
"graph");
154 APparams->set(
"compute global constants: temporaries",APparams->get(
"compute global constants: temporaries",
false));
155 APparams->set(
"compute global constants",APparams->get(
"compute global constants",
false));
157 const SC
dampingFactor = as<SC>(pL.get<
double>(
"sa: damping factor"));
158 const LO maxEigenIterations = as<LO>(pL.get<
int> (
"sa: eigenvalue estimate num iterations"));
159 const bool estimateMaxEigen = pL.get<
bool> (
"sa: calculate eigenvalue estimate");
160 const bool useAbsValueRowSum= pL.get<
bool> (
"sa: use rowsumabs diagonal scaling");
161 const bool doQRStep = pL.get<
bool>(
"tentative: calculate qr");
162 const bool enforceConstraints= pL.get<
bool>(
"sa: enforce constraints");
163 const SC userDefinedMaxEigen = as<SC>(pL.get<
double>(
"sa: max eigenvalue"));
164 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
165 double dTol = pL.get<
double>(
"sa: rowsumabs diagonal replacement tolerance");
166 const Magnitude diagonalReplacementTolerance = (dTol == as<double>(-1) ? Teuchos::ScalarTraits<Scalar>::eps()*100 : as<Magnitude>(pL.get<
double>(
"sa: rowsumabs diagonal replacement tolerance")));
167 const SC diagonalReplacementValue = as<SC>(pL.get<
double>(
"sa: rowsumabs diagonal replacement value"));
171 "MueLu::TentativePFactory::MakeTentative: cannot use 'enforce constraints' and 'calculate qr' at the same time");
176 Teuchos::RCP<Vector> invDiag;
177 if (userDefinedMaxEigen == -1.)
180 lambdaMax = A->GetMaxEigenvalueEstimate();
181 if (lambdaMax == -Teuchos::ScalarTraits<SC>::one() || estimateMaxEigen) {
182 GetOStream(
Statistics1) <<
"Calculating max eigenvalue estimate now (max iters = "<< maxEigenIterations <<
183 ( (useAbsValueRowSum) ?
", use rowSumAbs diagonal)" :
", use point diagonal)") << std::endl;
184 Coordinate stopTol = 1e-4;
185 if (useAbsValueRowSum) {
186 const bool returnReciprocal=
true;
187 const bool replaceSingleEntryRowWithZero=
true;
189 diagonalReplacementTolerance,
190 diagonalReplacementValue,
191 replaceSingleEntryRowWithZero);
193 "SaPFactory: eigenvalue estimate: diagonal reciprocal is null.");
197 A->SetMaxEigenvalueEstimate(lambdaMax);
199 GetOStream(
Statistics1) <<
"Using cached max eigenvalue estimate" << std::endl;
203 lambdaMax = userDefinedMaxEigen;
204 A->SetMaxEigenvalueEstimate(lambdaMax);
210 if (!useAbsValueRowSum)
212 else if (invDiag == Teuchos::null) {
213 GetOStream(
Runtime0) <<
"Using rowsumabs diagonal" << std::endl;
214 const bool returnReciprocal=
true;
215 const bool replaceSingleEntryRowWithZero=
true;
217 diagonalReplacementTolerance,
218 diagonalReplacementValue,
219 replaceSingleEntryRowWithZero);
220 TEUCHOS_TEST_FOR_EXCEPTION(invDiag.is_null(),
Exceptions::RuntimeError,
"SaPFactory: diagonal reciprocal is null.");
224 TEUCHOS_TEST_FOR_EXCEPTION(!std::isfinite(Teuchos::ScalarTraits<SC>::magnitude(omega)),
Exceptions::RuntimeError,
"Prolongator damping factor needs to be finite.");
227 finalP = Xpetra::IteratorOps<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP,
228 GetOStream(
Statistics2), std::string(
"MueLu::SaP-")+levelIDs, APparams);
229 if (enforceConstraints) SatisfyPConstraints( A, finalP);
237 if (!restrictionMode_) {
239 if(!finalP.is_null()) {std::ostringstream oss; oss <<
"P_" << coarseLevel.
GetLevelID(); finalP->setObjectLabel(oss.str());}
240 Set(coarseLevel,
"P", finalP);
242 APparams->set(
"graph", finalP);
243 Set(coarseLevel,
"AP reuse data", APparams);
246 if (Ptent->IsView(
"stridedMaps"))
247 finalP->CreateView(
"stridedMaps", Ptent);
250 RCP<ParameterList> params = rcp(
new ParameterList());
251 params->set(
"printLoadBalancingInfo",
true);
252 params->set(
"printCommInfo",
true);
262 if(!R.is_null()) {std::ostringstream oss; oss <<
"R_" << coarseLevel.
GetLevelID(); R->setObjectLabel(oss.str());}
265 Set(coarseLevel,
"R", R);
268 if (Ptent->IsView(
"stridedMaps"))
269 R->CreateView(
"stridedMaps", Ptent,
true);
272 RCP<ParameterList> params = rcp(
new ParameterList());
273 params->set(
"printLoadBalancingInfo",
true);
274 params->set(
"printCommInfo",
true);
297 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
300 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
301 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
302 LO nPDEs = A->GetFixedBlockSize();
303 Teuchos::ArrayRCP<Scalar> ConstraintViolationSum(nPDEs);
304 Teuchos::ArrayRCP<Scalar> Rsum(nPDEs);
305 Teuchos::ArrayRCP<size_t> nPositive(nPDEs);
306 for (
size_t k=0; k < (size_t) nPDEs; k++) ConstraintViolationSum[k] = zero;
307 for (
size_t k=0; k < (size_t) nPDEs; k++) Rsum[k] = zero;
308 for (
size_t k=0; k < (size_t) nPDEs; k++) nPositive[k] = 0;
311 for (
size_t i = 0; i < as<size_t>(P->getRowMap()->getNodeNumElements()); i++) {
313 Teuchos::ArrayView<const LocalOrdinal> indices;
314 Teuchos::ArrayView<const Scalar> vals1;
315 Teuchos::ArrayView< Scalar> vals;
316 P->getLocalRowView((LO) i, indices, vals1);
317 size_t nnz = indices.size();
318 if (nnz == 0)
continue;
320 vals = ArrayView<Scalar>(
const_cast<SC*
>(vals1.getRawPtr()), nnz);
323 bool checkRow =
true;
329 for (LO j = 0; j < indices.size(); j++) {
330 Rsum[ j%nPDEs ] += vals[j];
331 if (Teuchos::ScalarTraits<SC>::real(vals[j]) < Teuchos::ScalarTraits<SC>::real(zero)) {
332 ConstraintViolationSum[ j%nPDEs ] += vals[j];
336 if (Teuchos::ScalarTraits<SC>::real(vals[j]) != Teuchos::ScalarTraits<SC>::real(zero))
337 (nPositive[ j%nPDEs])++;
339 if (Teuchos::ScalarTraits<SC>::real(vals[j]) > Teuchos::ScalarTraits<SC>::real(1.00001 )) {
340 ConstraintViolationSum[ j%nPDEs ] += (vals[j] - one);
350 for (
size_t k=0; k < (size_t) nPDEs; k++) {
352 if (Teuchos::ScalarTraits<SC>::real(Rsum[ k ]) < Teuchos::ScalarTraits<SC>::magnitude(zero)) {
353 ConstraintViolationSum[k] += (-Rsum[k]);
355 else if (Teuchos::ScalarTraits<SC>::real(Rsum[ k ]) > Teuchos::ScalarTraits<SC>::magnitude(1.00001)) {
356 ConstraintViolationSum[k] += (one - Rsum[k]);
361 for (
size_t k=0; k < (size_t) nPDEs; k++) {
362 if (Teuchos::ScalarTraits<SC>::magnitude(ConstraintViolationSum[ k ]) != Teuchos::ScalarTraits<SC>::magnitude(zero))
367 for (LO j = 0; j < indices.size(); j++) {
368 if (Teuchos::ScalarTraits<SC>::real(vals[j]) > Teuchos::ScalarTraits<SC>::real(zero)) {
369 vals[j] += (ConstraintViolationSum[j%nPDEs]/ (as<Scalar>(nPositive[j%nPDEs])));
372 for (
size_t k=0; k < (size_t) nPDEs; k++) ConstraintViolationSum[k] = zero;
374 for (
size_t k=0; k < (size_t) nPDEs; k++) Rsum[k] = zero;
375 for (
size_t k=0; k < (size_t) nPDEs; k++) nPositive[k] = 0;
382 #endif // MUELU_SAPFACTORY_DEF_HPP void SatisfyPConstraints(RCP< Matrix > A, RCP< Matrix > &P) const
Enforce constraints on prolongator.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
scalar_type dampingFactor
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetLumpedMatrixDiagonal(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Print even more statistics.
int GetLevelID() const
Return level number.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MueLu::DefaultScalar Scalar
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
#define SET_VALID_ENTRY(name)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Exception throws to report errors in the internal logical of the program.
static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixDiagonalInverse(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
static Scalar PowerMethod(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)