46#ifndef MUELU_REITZINGERPFACTORY_DEF_HPP
47#define MUELU_REITZINGERPFACTORY_DEF_HPP
49#include <Xpetra_MapFactory.hpp>
50#include <Xpetra_Map.hpp>
51#include <Xpetra_CrsMatrix.hpp>
52#include <Xpetra_Matrix.hpp>
53#include <Xpetra_MatrixMatrix.hpp>
54#include <Xpetra_MultiVector.hpp>
55#include <Xpetra_MultiVectorFactory.hpp>
56#include <Xpetra_VectorFactory.hpp>
57#include <Xpetra_Import.hpp>
58#include <Xpetra_ImportUtils.hpp>
59#include <Xpetra_ImportFactory.hpp>
60#include <Xpetra_CrsMatrixWrap.hpp>
61#include <Xpetra_StridedMap.hpp>
62#include <Xpetra_StridedMapFactory.hpp>
63#include <Xpetra_IO.hpp>
67#include "MueLu_Aggregates.hpp"
68#include "MueLu_AmalgamationFactory.hpp"
69#include "MueLu_AmalgamationInfo.hpp"
70#include "MueLu_CoarseMapFactory.hpp"
73#include "MueLu_NullspaceFactory.hpp"
74#include "MueLu_PerfUtils.hpp"
75#include "MueLu_Utilities.hpp"
79 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 RCP<ParameterList> validParamList = rcp(
new ParameterList());
83#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
90 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
91 validParamList->set< RCP<const FactoryBase> >(
"D0", Teuchos::null,
"Generating factory of the matrix D0");
92 validParamList->set< RCP<const FactoryBase> >(
"NodeAggMatrix", Teuchos::null,
"Generating factory of the matrix NodeAggMatrix");
93 validParamList->set< RCP<const FactoryBase> >(
"Pnodal", Teuchos::null,
"Generating factory of the matrix P");
94 validParamList->set< RCP<const FactoryBase> >(
"NodeImporter", Teuchos::null,
"Generating factory of the matrix NodeImporter");
97 ParameterList norecurse;
98 norecurse.disableRecursiveValidation();
99 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
101 return validParamList;
104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
106 Input(fineLevel,
"A");
107 Input(fineLevel,
"D0");
108 Input(fineLevel,
"NodeAggMatrix");
109 Input(coarseLevel,
"NodeAggMatrix");
110 Input(coarseLevel,
"Pnodal");
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 return BuildP(fineLevel, coarseLevel);
120 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
123 using Teuchos::arcp_const_cast;
124 using MT =
typename Teuchos::ScalarTraits<SC>::magnitudeType;
125 using XMM = Xpetra::MatrixMatrix<SC,LO,GO,NO>;
126 Teuchos::FancyOStream& out0=GetBlackHole();
127 const ParameterList& pL = GetParameterList();
129 bool update_communicators = pL.get<
bool>(
"repartition: enable") && pL.get<
bool>(
"repartition: use subcommunicators");
132 bool nodal_p_is_all_ones = !pL.get<
bool>(
"tentative: constant column sums") && !pL.get<
bool>(
"tentative: calculate qr");
134 RCP<Matrix> EdgeMatrix = Get< RCP<Matrix> > (fineLevel,
"A");
135 RCP<Matrix> D0 = Get< RCP<Matrix> > (fineLevel,
"D0");
136 RCP<Matrix> NodeMatrix = Get< RCP<Matrix> > (fineLevel,
"NodeAggMatrix");
137 RCP<Matrix> Pn = Get< RCP<Matrix> > (coarseLevel,
"Pnodal");
139 const GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
140 const LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
143 RCP<Operator> CoarseNodeMatrix = Get< RCP<Operator> >(coarseLevel,
"NodeAggMatrix");
144 int MyPID = EdgeMatrix.is_null()? -1 : EdgeMatrix->getRowMap()->getComm()->getRank();
147 RCP<ParameterList> mm_params = rcp(
new ParameterList);;
148 if(pL.isSublist(
"matrixmatrix: kernel params"))
149 mm_params->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
153 if(!nodal_p_is_all_ones) {
155 GetOStream(
Runtime0) <<
"ReitzingerPFactory::BuildP(): Assuming Pn is not normalized" << std::endl;
156 RCP<Matrix> Pn_old = Pn;
158 Pn = Xpetra::MatrixFactory<SC,LO,GO,NO>::Build(Pn->getCrsGraph());
159 Pn->setAllToScalar(Teuchos::ScalarTraits<SC>::one());
160 Pn->fillComplete(Pn->getDomainMap(),Pn->getRangeMap());
164 GetOStream(
Runtime0) <<
"ReitzingerPFactory::BuildP(): Assuming Pn is normalized" << std::endl;
172 RCP<Matrix> D0_Pn, PnT_D0T, D0_Pn_nonghosted;
173 Teuchos::Array<int> D0_Pn_col_pids;
177 D0_Pn = XMM::Multiply(*D0,
false,*Pn,
false,dummy,out0,
true,
true,
"D0*Pn",mm_params);
180 if(!mm_params.is_null()) mm_params->remove(
"importer",
false);
183 D0_Pn_nonghosted = D0_Pn;
186 if(!D0_Pn->getCrsGraph()->getImporter().is_null()) {
187 Xpetra::ImportUtils<LO,GO,NO> utils;
188 utils.getPids(*D0_Pn->getCrsGraph()->getImporter(),D0_Pn_col_pids,
false);
191 D0_Pn_col_pids.resize(D0_Pn->getCrsGraph()->getColMap()->getLocalNumElements(),MyPID);
208 if(!PnT_D0T->getCrsGraph()->getImporter().is_null()) {
209 RCP<const Import> Importer = PnT_D0T->getCrsGraph()->getImporter();
210 RCP<const CrsMatrix> D0_Pn_crs = rcp_dynamic_cast<const CrsMatrixWrap>(D0_Pn)->getCrsMatrix();
211 RCP<Matrix> D0_Pn_new = rcp(
new CrsMatrixWrap(CrsMatrixFactory::Build(D0_Pn_crs,*Importer,D0_Pn->getDomainMap(),Importer->getTargetMap())));
214 if(!D0_Pn->getCrsGraph()->getImporter().is_null()) {
215 Xpetra::ImportUtils<LO,GO,NO> utils;
216 utils.getPids(*D0_Pn->getCrsGraph()->getImporter(),D0_Pn_col_pids,
false);
219 D0_Pn_col_pids.resize(D0_Pn->getCrsGraph()->getColMap()->getLocalNumElements(),MyPID);
225 ArrayView<const LO> colind_E, colind_N;
226 ArrayView<const SC> values_E, values_N;
228 size_t Ne=EdgeMatrix->getLocalNumRows();
229 size_t Nn=NodeMatrix->getLocalNumRows();
232 size_t max_edges = (NodeMatrix->getLocalNumEntries() + Nn +1) / 2;
233 ArrayRCP<size_t> D0_rowptr(Ne+1);
234 ArrayRCP<LO> D0_colind(max_edges);
235 ArrayRCP<SC> D0_values(max_edges);
239 LO Nnc = PnT_D0T->getRowMap()->getLocalNumElements();
241 for(LO i=0; i<(LO)Nnc; i++) {
245 using value_type = bool;
246 std::map<LO, value_type> ce_map;
249 PnT_D0T->getLocalRowView(i,colind_E,values_E);
251 for(LO j=0; j<(LO)colind_E.size(); j++) {
262 GO edge_gid = PnT_D0T->getColMap()->getGlobalElement(colind_E[j]);
263 LO j_row = D0_Pn->getRowMap()->getLocalElement(edge_gid);
265 D0_Pn->getLocalRowView(j_row,colind_N,values_N);
268 if(colind_N.size() != 2)
continue;
270 pid0 = D0_Pn_col_pids[colind_N[0]];
271 pid1 = D0_Pn_col_pids[colind_N[1]];
277 bool zero_matches = pid0 == MyPID;
278 bool one_matches = pid1 == MyPID;
279 bool keep_shared_edge =
false, own_both_nodes =
false;
280 if(zero_matches && one_matches) {own_both_nodes=
true;}
282 int sum_is_odd = (pid0 + pid1) % 2;
283 int i_am_smaller = MyPID == std::min(pid0,pid1);
284 if(sum_is_odd && i_am_smaller) keep_shared_edge=
true;
285 if(!sum_is_odd && !i_am_smaller) keep_shared_edge=
true;
288 if(!keep_shared_edge && !own_both_nodes)
continue;
295 for(LO k=0; k<(LO)colind_N.size(); k++) {
296 LO my_colind = colind_N[k];
297 if(my_colind!=LO_INVALID && ((keep_shared_edge && my_colind != i) || (own_both_nodes && my_colind > i)) ) {
298 ce_map.emplace(std::make_pair(my_colind,
true));
305 for(
auto iter=ce_map.begin(); iter != ce_map.end(); iter++) {
306 LO col = iter->first;
313 D0_colind[current] = i;
314 D0_values[current] = -1;
316 D0_colind[current] = col;
317 D0_values[current] = 1;
319 D0_rowptr[current / 2] = current;
324 LO num_coarse_edges = current / 2;
325 D0_rowptr.resize(num_coarse_edges+1);
326 D0_colind.resize(current);
327 D0_values.resize(current);
331 TEUCHOS_TEST_FOR_EXCEPTION( (num_coarse_edges > 0 && CoarseNodeMatrix.is_null()) ||
332 (num_coarse_edges == 0 && !CoarseNodeMatrix.is_null())
338 RCP<const Map> ownedCoarseEdgeMap = Xpetra::MapFactory<LO,GO,NO>::Build(EdgeMatrix->getRowMap()->lib(), GO_INVALID, num_coarse_edges,EdgeMatrix->getRowMap()->getIndexBase(),EdgeMatrix->getRowMap()->getComm());
342 RCP<const Map> ownedCoarseNodeMap = Pn->getDomainMap();
343 RCP<const Map> ownedPlusSharedCoarseNodeMap = D0_Pn->getCrsGraph()->getColMap();
346 RCP<CrsMatrix> D0_coarse;
351 D0_coarse = CrsMatrixFactory::Build(ownedCoarseEdgeMap,ownedPlusSharedCoarseNodeMap,0);
352 TEUCHOS_TEST_FOR_EXCEPTION(D0_coarse.is_null(),
Exceptions::RuntimeError,
"MueLu::ReitzingerPFactory: CrsMatrixFatory failed.");
358 D0_coarse->allocateAllValues(current, ia, ja, val);
359 std::copy(D0_rowptr.begin(),D0_rowptr.end(),ia.begin());
360 std::copy(D0_colind.begin(),D0_colind.end(),ja.begin());
361 std::copy(D0_values.begin(),D0_values.end(),val.begin());
362 D0_coarse->setAllValues(ia, ja, val);
367 printf(
"[%d] D0: ia.size() = %d ja.size() = %d\n",MyPID,(
int)ia.size(),(
int)ja.size());
368 printf(
"[%d] D0: ia :",MyPID);
369 for(
int i=0; i<(int)ia.size(); i++)
370 printf(
"%d ",(
int)ia[i]);
371 printf(
"\n[%d] D0: global ja :",MyPID);
372 for(
int i=0; i<(int)ja.size(); i++)
373 printf(
"%d ",(
int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i]));
374 printf(
"\n[%d] D0: local ja :",MyPID);
375 for(
int i=0; i<(int)ja.size(); i++)
376 printf(
"%d ",(
int)ja[i]);
379 sprintf(fname,
"D0_global_ja_%d_%d.dat",MyPID,fineLevel.
GetLevelID());
380 FILE * f = fopen(fname,
"w");
381 for(
int i=0; i<(int)ja.size(); i++)
382 fprintf(f,
"%d ",(
int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i]));
385 sprintf(fname,
"D0_local_ja_%d_%d.dat",MyPID,fineLevel.
GetLevelID());
386 f = fopen(fname,
"w");
387 for(
int i=0; i<(int)ja.size(); i++)
388 fprintf(f,
"%d ",(
int)ja[i]);
393 D0_coarse->expertStaticFillComplete(ownedCoarseNodeMap,ownedCoarseEdgeMap);
395 RCP<Matrix> D0_coarse_m = rcp(
new CrsMatrixWrap(D0_coarse));
396 RCP<Teuchos::FancyOStream> fout = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
410 RCP<Matrix> Pn_D0cT = XMM::Multiply(*Pn,
false,*D0_coarse_m,
true,dummy,out0,
true,
true,
"Pn*D0c'",mm_params);
413 if(!mm_params.is_null()) mm_params->remove(
"importer",
false);
415 Pe = XMM::Multiply(*D0,
false,*Pn_D0cT,
false,dummy,out0,
true,
true,
"D0*(Pn*D0c')",mm_params);
425 SC one = Teuchos::ScalarTraits<SC>::one();
426 MT two = 2*Teuchos::ScalarTraits<MT>::one();
427 SC zero = Teuchos::ScalarTraits<SC>::zero();
430 RCP<const CrsMatrix> Pe_crs = rcp_dynamic_cast<const CrsMatrixWrap>(Pe)->getCrsMatrix();
431 TEUCHOS_TEST_FOR_EXCEPTION(Pe_crs.is_null(),
Exceptions::RuntimeError,
"MueLu::ReitzingerPFactory: Pe is not a crs matrix.");
432 ArrayRCP<const size_t > rowptr_const;
433 ArrayRCP<const LO> colind_const;
434 ArrayRCP<const SC> values_const;
435 Pe_crs->getAllValues(rowptr_const,colind_const,values_const);
436 ArrayRCP<size_t> rowptr = arcp_const_cast<size_t>(rowptr_const);
437 ArrayRCP<LO> colind = arcp_const_cast<LO>(colind_const);
438 ArrayRCP<SC> values = arcp_const_cast<SC>(values_const);
440 LO lower = rowptr[0];
441 for(LO i=0; i<(LO)Ne; i++) {
442 for(
size_t j=lower; j<rowptr[i+1]; j++) {
443 if (values[j] == one || values[j] == neg_one || values[j] == zero) {
447 colind[ct] = colind[j];
448 values[ct] = values[j] / two;
458 rcp_const_cast<CrsMatrix>(Pe_crs)->setAllValues(rowptr,colind,values);
460 Pe->fillComplete(Pe->getDomainMap(),Pe->getRangeMap());
464 CheckCommutingProperty(*Pe,*D0_coarse_m,*D0,*Pn);
469 if(update_communicators) {
471 RCP<const Teuchos::Comm<int> > newComm;
472 if(!CoarseNodeMatrix.is_null()) newComm = CoarseNodeMatrix->getDomainMap()->getComm();
473 RCP<const Map> newMap = Xpetra::MapFactory<LO,GO,NO>::copyMapWithNewComm(D0_coarse_m->getRowMap(),newComm);
474 D0_coarse_m->removeEmptyProcessesInPlace(newMap);
477 if(newMap.is_null()) D0_coarse_m = Teuchos::null;
479 Set(coarseLevel,
"InPlaceMap",newMap);
483 Set(coarseLevel,
"P",Pe);
484 Set(coarseLevel,
"Ptent",Pe);
486 Set(coarseLevel,
"D0",D0_coarse_m);
493 int numProcs = Pe->getRowMap()->getComm()->getSize();
496 sprintf(fname,
"Pe_%d_%d.mat",numProcs,fineLevel.
GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*Pe);
497 sprintf(fname,
"Pn_%d_%d.mat",numProcs,fineLevel.
GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*Pn);
498 sprintf(fname,
"D0c_%d_%d.mat",numProcs,fineLevel.
GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*D0_coarse_m);
499 sprintf(fname,
"D0f_%d_%d.mat",numProcs,fineLevel.
GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*D0);
505 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
507 CheckCommutingProperty(
const Matrix & Pe,
const Matrix & D0_c,
const Matrix& D0_f,
const Matrix & Pn)
const {
509 using XMM = Xpetra::MatrixMatrix<SC,LO,GO,NO>;
510 using MT =
typename Teuchos::ScalarTraits<SC>::magnitudeType;
511 SC one = Teuchos::ScalarTraits<SC>::one();
512 SC zero = Teuchos::ScalarTraits<SC>::zero();
515 Teuchos::FancyOStream &out0=GetBlackHole();
516 RCP<Matrix> left = XMM::Multiply(Pe,
false,D0_c,
false,dummy,out0);
517 RCP<Matrix> right = XMM::Multiply(D0_f,
false,Pn,
false,dummy,out0);
520 RCP<CrsMatrix> sum_c = CrsMatrixFactory::Build(left->getRowMap(),left->getLocalMaxNumRowEntries()+right->getLocalMaxNumRowEntries());
521 RCP<Matrix> summation = rcp(
new CrsMatrixWrap(sum_c));
522 XMM::TwoMatrixAdd(*left,
false, one, *summation, zero);
523 XMM::TwoMatrixAdd(*right,
false, -one, *summation, one);
525 MT norm = summation->getFrobeniusNorm();
526 GetOStream(
Statistics0) <<
"CheckCommutingProperty: ||Pe D0_c - D0_f Pn || = "<<norm<<std::endl;
538#define MUELU_REITZINGERPFACTORY_SHORT
#define SET_VALID_ENTRY(name)
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
int GetLevelID() const
Return level number.
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const NoFactory * get()
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void CheckCommutingProperty(const Matrix &Pe, const Matrix &D0_c, const Matrix &D0_f, const Matrix &Pn) const
Utility method.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
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)
Namespace for MueLu classes and methods.
@ Final
Keep data only for this run. Used to keep data useful for Hierarchy::Iterate(). Data will be deleted ...
@ UserData
User data are always kept. This flag is set automatically when Level::Set("data", data) is used....
@ Runtime0
One-liner description of what is happening.
@ Statistics0
Print statistics that do not involve significant additional computation.