46 #ifndef MUELU_REITZINGERPFACTORY_DEF_HPP
47 #define MUELU_REITZINGERPFACTORY_DEF_HPP
49 #include <Xpetra_MapFactory.hpp>
50 #include <Xpetra_Map.hpp>
55 #include <Xpetra_MultiVectorFactory.hpp>
58 #include <Xpetra_ImportUtils.hpp>
60 #include <Xpetra_CrsMatrixWrap.hpp>
61 #include <Xpetra_StridedMap.hpp>
62 #include <Xpetra_StridedMapFactory.hpp>
69 #include "MueLu_Utilities.hpp"
73 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
77 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
82 #undef SET_VALID_ENTRY
86 validParamList->
set<
RCP<const FactoryBase> >(
"NodeAggMatrix", Teuchos::null,
"Generating factory of the matrix NodeAggMatrix");
88 validParamList->
set<
RCP<const FactoryBase> >(
"NodeImporter", Teuchos::null,
"Generating factory of the matrix NodeImporter");
92 norecurse.disableRecursiveValidation();
93 validParamList->
set<
ParameterList>(
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
95 return validParamList;
98 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 Input(fineLevel,
"A");
101 Input(fineLevel,
"D0");
102 Input(fineLevel,
"NodeAggMatrix");
103 Input(coarseLevel,
"NodeAggMatrix");
104 Input(coarseLevel,
"Pnodal");
108 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
110 return BuildP(fineLevel, coarseLevel);
113 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
116 using Teuchos::arcp_const_cast;
122 bool update_communicators = pL.
get<
bool>(
"repartition: enable") && pL.
get<
bool>(
"repartition: use subcommunicators");
125 bool nodal_p_is_all_ones = !pL.
get<
bool>(
"tentative: constant column sums") && !pL.
get<
bool>(
"tentative: calculate qr");
127 RCP<Matrix> EdgeMatrix = Get<RCP<Matrix> >(fineLevel,
"A");
128 RCP<Matrix> D0 = Get<RCP<Matrix> >(fineLevel,
"D0");
129 RCP<Matrix> NodeMatrix = Get<RCP<Matrix> >(fineLevel,
"NodeAggMatrix");
130 RCP<Matrix> Pn = Get<RCP<Matrix> >(coarseLevel,
"Pnodal");
136 RCP<Operator> CoarseNodeMatrix = Get<RCP<Operator> >(coarseLevel,
"NodeAggMatrix");
137 int MyPID = EdgeMatrix.
is_null() ? -1 : EdgeMatrix->getRowMap()->getComm()->getRank();
142 if (pL.
isSublist(
"matrixmatrix: kernel params"))
143 mm_params->
sublist(
"matrixmatrix: kernel params") = pL.
sublist(
"matrixmatrix: kernel params");
146 if (!nodal_p_is_all_ones) {
148 GetOStream(
Runtime0) <<
"ReitzingerPFactory::BuildP(): Assuming Pn is not normalized" << std::endl;
153 Pn->fillComplete(Pn->getDomainMap(), Pn->getRangeMap());
156 GetOStream(
Runtime0) <<
"ReitzingerPFactory::BuildP(): Assuming Pn is normalized" << std::endl;
169 D0_Pn = XMM::Multiply(*D0,
false, *Pn,
false, dummy, out0,
true,
true,
"D0*Pn", mm_params);
172 if (!mm_params.
is_null()) mm_params->
remove(
"importer",
false);
175 D0_Pn_nonghosted = D0_Pn;
178 if (!D0_Pn->getCrsGraph()->getImporter().
is_null()) {
179 Xpetra::ImportUtils<LO, GO, NO> utils;
180 utils.getPids(*D0_Pn->getCrsGraph()->getImporter(), D0_Pn_col_pids,
false);
182 D0_Pn_col_pids.
resize(D0_Pn->getCrsGraph()->getColMap()->getLocalNumElements(), MyPID);
197 if (!PnT_D0T->getCrsGraph()->getImporter().
is_null()) {
200 RCP<Matrix> D0_Pn_new =
rcp(
new CrsMatrixWrap(CrsMatrixFactory::Build(D0_Pn_crs, *Importer, D0_Pn->getDomainMap(), Importer->getTargetMap())));
203 if (!D0_Pn->getCrsGraph()->getImporter().is_null()) {
204 Xpetra::ImportUtils<LO, GO, NO> utils;
205 utils.getPids(*D0_Pn->getCrsGraph()->getImporter(), D0_Pn_col_pids,
false);
207 D0_Pn_col_pids.
resize(D0_Pn->getCrsGraph()->getColMap()->getLocalNumElements(), MyPID);
215 size_t Ne = EdgeMatrix->getLocalNumRows();
216 size_t Nn = NodeMatrix->getLocalNumRows();
219 size_t max_edges = (NodeMatrix->getLocalNumEntries() + Nn + 1) / 2;
226 LO Nnc = PnT_D0T->getRowMap()->getLocalNumElements();
230 RCP<const Map> ownedPlusSharedCoarseNodeMap = D0_Pn->getCrsGraph()->getColMap();
232 for (
LO i = 0; i < (
LO)Nnc; i++) {
233 LO local_column_i = ownedPlusSharedCoarseNodeMap->getLocalElement(PnT_D0T->getRowMap()->getGlobalElement(i));
236 using value_type = bool;
237 std::map<LO, value_type> ce_map;
240 PnT_D0T->getLocalRowView(i, colind_E, values_E);
242 for (
LO j = 0; j < (
LO)colind_E.
size(); j++) {
252 GO edge_gid = PnT_D0T->getColMap()->getGlobalElement(colind_E[j]);
253 LO j_row = D0_Pn->getRowMap()->getLocalElement(edge_gid);
255 D0_Pn->getLocalRowView(j_row, colind_N, values_N);
258 if (colind_N.
size() != 2)
continue;
260 pid0 = D0_Pn_col_pids[colind_N[0]];
261 pid1 = D0_Pn_col_pids[colind_N[1]];
267 bool zero_matches = pid0 == MyPID;
268 bool one_matches = pid1 == MyPID;
269 bool keep_shared_edge =
false, own_both_nodes =
false;
270 if (zero_matches && one_matches) {
271 own_both_nodes =
true;
273 bool sum_is_even = (pid0 + pid1) % 2 == 0;
274 bool i_am_smaller = MyPID == std::min(pid0, pid1);
275 if (sum_is_even && i_am_smaller) keep_shared_edge =
true;
276 if (!sum_is_even && !i_am_smaller) keep_shared_edge =
true;
279 if (!keep_shared_edge && !own_both_nodes)
continue;
285 for (
LO k = 0; k < (
LO)colind_N.size(); k++) {
286 LO my_colind = colind_N[k];
287 if (my_colind != LO_INVALID && ((keep_shared_edge && my_colind != local_column_i) || (own_both_nodes && my_colind > local_column_i))) {
288 ce_map.emplace(std::make_pair(my_colind,
true));
294 for (
auto iter = ce_map.begin();
iter != ce_map.end();
iter++) {
296 if (col == local_column_i) {
301 D0_colind[current] = local_column_i;
302 D0_values[current] = -1;
304 D0_colind[current] = col;
305 D0_values[current] = 1;
307 D0_rowptr[current / 2] = current;
312 LO num_coarse_edges = current / 2;
313 D0_rowptr.
resize(num_coarse_edges + 1);
314 D0_colind.
resize(current);
315 D0_values.
resize(current);
327 D0_coarse = CrsMatrixFactory::Build(ownedCoarseEdgeMap, ownedPlusSharedCoarseNodeMap, 0);
334 D0_coarse->allocateAllValues(current, ia, ja, val);
335 std::copy(D0_rowptr.
begin(), D0_rowptr.
end(), ia.begin());
338 D0_coarse->setAllValues(ia, ja, val);
343 printf(
"[%d] D0: ia.size() = %d ja.size() = %d\n",MyPID,(
int)ia.size(),(int)ja.
size());
344 printf(
"[%d] D0: ia :",MyPID);
345 for(
int i=0; i<(int)ia.size(); i++)
346 printf(
"%d ",(
int)ia[i]);
347 printf(
"\n[%d] D0: global ja :",MyPID);
348 for(
int i=0; i<(int)ja.
size(); i++)
349 printf(
"%d ",(
int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i]));
350 printf(
"\n[%d] D0: local ja :",MyPID);
351 for(
int i=0; i<(int)ja.
size(); i++)
352 printf(
"%d ",(
int)ja[i]);
355 sprintf(fname,
"D0_global_ja_%d_%d.dat",MyPID,fineLevel.GetLevelID());
356 FILE * f = fopen(fname,
"w");
357 for(
int i=0; i<(int)ja.
size(); i++)
358 fprintf(f,
"%d ",(
int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i]));
361 sprintf(fname,
"D0_local_ja_%d_%d.dat",MyPID,fineLevel.GetLevelID());
362 f = fopen(fname,
"w");
363 for(
int i=0; i<(int)ja.
size(); i++)
364 fprintf(f,
"%d ",(
int)ja[i]);
369 D0_coarse->expertStaticFillComplete(ownedCoarseNodeMap, ownedCoarseEdgeMap);
385 RCP<Matrix> Pn_D0cT = XMM::Multiply(*Pn,
false, *D0_coarse_m,
true, dummy, out0,
true,
true,
"Pn*D0c'", mm_params);
388 if (!mm_params.
is_null()) mm_params->
remove(
"importer",
false);
390 Pe = XMM::Multiply(*D0,
false, *Pn_D0cT,
false, dummy, out0,
true,
true,
"D0*(Pn*D0c')", mm_params);
410 Pe_crs->getAllValues(rowptr_const, colind_const, values_const);
415 LO lower = rowptr[0];
416 for (
LO i = 0; i < (
LO)Ne; i++) {
417 for (
size_t j = lower; j < rowptr[i + 1]; j++) {
418 if (values[j] == one || values[j] == neg_one || values[j] == zero) {
421 colind[ct] = colind[j];
422 values[ct] = values[j] / two;
426 lower = rowptr[i + 1];
432 rcp_const_cast<CrsMatrix>(Pe_crs)->setAllValues(rowptr, colind, values);
434 Pe->fillComplete(Pe->getDomainMap(), Pe->getRangeMap());
438 CheckCommutingProperty(*Pe, *D0_coarse_m, *D0, *Pn);
443 if (update_communicators) {
446 if (!CoarseNodeMatrix.
is_null()) newComm = CoarseNodeMatrix->getDomainMap()->getComm();
448 D0_coarse_m->removeEmptyProcessesInPlace(newMap);
451 if (newMap.
is_null()) D0_coarse_m = Teuchos::null;
453 Set(coarseLevel,
"InPlaceMap", newMap);
457 Set(coarseLevel,
"P", Pe);
458 Set(coarseLevel,
"Ptent", Pe);
460 Set(coarseLevel,
"D0", D0_coarse_m);
467 int numProcs = Pe->getRowMap()->getComm()->getSize();
479 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
490 RCP<Matrix> left = XMM::Multiply(Pe,
false, D0_c,
false, dummy, out0);
491 RCP<Matrix> right = XMM::Multiply(D0_f,
false, Pn,
false, dummy, out0);
494 RCP<CrsMatrix> sum_c = CrsMatrixFactory::Build(left->getRowMap(), left->getLocalMaxNumRowEntries() + right->getLocalMaxNumRowEntries());
496 XMM::TwoMatrixAdd(*left,
false, one, *summation, zero);
497 XMM::TwoMatrixAdd(*right,
false, -one, *summation, one);
499 MT norm = summation->getFrobeniusNorm();
500 GetOStream(
Statistics0) <<
"CheckCommutingProperty: ||Pe D0_c - D0_f Pn || = " << norm << std::endl;
507 #define MUELU_REITZINGERPFACTORY_SHORT
508 #endif // MUELU_REITZINGERPFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
One-liner description of what is happening.
static const NoFactory * get()
void resize(const size_type n, const T &val=T())
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Print statistics that do not involve significant additional computation.
bool remove(std::string const &name, bool throwIfNotExists=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
bool isSublist(const std::string &name) const
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Keep data only for this run. Used to keep data useful for Hierarchy::Iterate(). Data will be deleted ...
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 resize(size_type new_size, const value_type &x=value_type())
void CheckCommutingProperty(const Matrix &Pe, const Matrix &D0_c, const Matrix &D0_f, const Matrix &Pn) const
Utility method.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Exception throws to report errors in the internal logical of the program.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.