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)) 86 #undef SET_VALID_ENTRY 88 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
89 validParamList->set< RCP<const FactoryBase> >(
"D0", Teuchos::null,
"Generating factory of the matrix D0");
90 validParamList->set< RCP<const FactoryBase> >(
"NodeMatrix", Teuchos::null,
"Generating factory of the matrix NodeMatrix");
91 validParamList->set< RCP<const FactoryBase> >(
"Pnodal", Teuchos::null,
"Generating factory of the matrix P");
94 ParameterList norecurse;
95 norecurse.disableRecursiveValidation();
96 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
98 return validParamList;
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
103 Input(fineLevel,
"A");
104 Input(fineLevel,
"D0");
105 Input(fineLevel,
"NodeMatrix");
106 Input(coarseLevel,
"Pnodal");
110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
112 return BuildP(fineLevel, coarseLevel);
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
118 using MT =
typename Teuchos::ScalarTraits<SC>::magnitudeType;
119 using XMM = Xpetra::MatrixMatrix<SC,LO,GO,NO>;
120 Teuchos::FancyOStream& out0=GetBlackHole();
121 const ParameterList& pL = GetParameterList();
123 RCP<Matrix> EdgeMatrix = Get< RCP<Matrix> > (fineLevel,
"A");
124 RCP<Matrix> D0 = Get< RCP<Matrix> > (fineLevel,
"D0");
125 RCP<Matrix> NodeMatrix = Get< RCP<Matrix> > (fineLevel,
"NodeMatrix");
126 RCP<Matrix> Pn = Get< RCP<Matrix> > (coarseLevel,
"Pnodal");
127 const GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
128 const LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
129 int MyPID = D0->getRowMap()->getComm()->getRank();
132 RCP<ParameterList> mm_params = rcp(
new ParameterList);;
133 if(pL.isSublist(
"matrixmatrix: kernel params"))
134 mm_params->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
142 RCP<Matrix> D0_Pn, PnT_D0T, D0_Pn_nonghosted;
143 Teuchos::Array<int> D0_Pn_col_pids;
147 D0_Pn = XMM::Multiply(*D0,
false,*Pn,
false,dummy,out0,
true,
true,
"D0*Pn",mm_params);
150 D0_Pn_nonghosted = D0_Pn;
153 if(!D0_Pn->getCrsGraph()->getImporter().is_null()) {
154 Xpetra::ImportUtils<LO,GO,NO> utils;
155 utils.getPids(*D0_Pn->getCrsGraph()->getImporter(),D0_Pn_col_pids,
false);
158 D0_Pn_col_pids.resize(D0_Pn->getCrsGraph()->getColMap()->getNodeNumElements(),MyPID);
174 if(!PnT_D0T->getCrsGraph()->getImporter().is_null()) {
175 RCP<const Import> Importer = PnT_D0T->getCrsGraph()->getImporter();
176 RCP<const CrsMatrix> D0_Pn_crs = rcp_dynamic_cast<
const CrsMatrixWrap>(D0_Pn)->getCrsMatrix();
177 RCP<Matrix> D0_Pn_new = rcp(
new CrsMatrixWrap(CrsMatrixFactory::Build(D0_Pn_crs,*Importer,D0_Pn->getDomainMap(),Importer->getTargetMap())));
180 if(!D0_Pn->getCrsGraph()->getImporter().is_null()) {
181 Xpetra::ImportUtils<LO,GO,NO> utils;
182 utils.getPids(*D0_Pn->getCrsGraph()->getImporter(),D0_Pn_col_pids,
false);
185 D0_Pn_col_pids.resize(D0_Pn->getCrsGraph()->getColMap()->getNodeNumElements(),MyPID);
191 ArrayView<const LO> colind_E, colind_N;
192 ArrayView<const SC> values_E, values_N;
194 size_t Ne=EdgeMatrix->getNodeNumRows();
195 size_t Nn=NodeMatrix->getNodeNumRows();
198 size_t max_edges = (NodeMatrix->getNodeNumEntries() + Nn +1) / 2;
199 ArrayRCP<size_t> D0_rowptr(Ne+1);
200 ArrayRCP<LO> D0_colind(max_edges);
201 ArrayRCP<SC> D0_values(max_edges);
205 LO Nnc = PnT_D0T->getRowMap()->getNodeNumElements();
207 for(LO i=0; i<(LO)Nnc; i++) {
211 using value_type = bool;
212 std::map<LO, value_type> ce_map;
215 PnT_D0T->getLocalRowView(i,colind_E,values_E);
217 for(LO j=0; j<(LO)colind_E.size(); j++) {
228 GO edge_gid = PnT_D0T->getColMap()->getGlobalElement(colind_E[j]);
229 LO j_row = D0_Pn->getRowMap()->getLocalElement(edge_gid);
231 D0_Pn->getLocalRowView(j_row,colind_N,values_N);
234 if(colind_N.size() != 2)
continue;
236 pid0 = D0_Pn_col_pids[colind_N[0]];
237 pid1 = D0_Pn_col_pids[colind_N[1]];
243 bool zero_matches = pid0 == MyPID;
244 bool one_matches = pid1 == MyPID;
245 bool keep_shared_edge =
false, own_both_nodes =
false;
246 if(zero_matches && one_matches) {own_both_nodes=
true;}
248 int sum_is_odd = (pid0 + pid1) % 2;
249 int i_am_smaller = MyPID == std::min(pid0,pid1);
250 if(sum_is_odd && i_am_smaller) keep_shared_edge=
true;
251 if(!sum_is_odd && !i_am_smaller) keep_shared_edge=
true;
254 if(!keep_shared_edge && !own_both_nodes)
continue;
261 for(LO k=0; k<(LO)colind_N.size(); k++) {
262 LO my_colind = colind_N[k];
263 if(my_colind!=LO_INVALID && ((keep_shared_edge && my_colind != i) || (own_both_nodes && my_colind > i)) ) {
264 ce_map.emplace(std::make_pair(my_colind,
true));
271 for(
auto iter=ce_map.begin(); iter != ce_map.end(); iter++) {
272 LO col = iter->first;
279 D0_colind[current] = i;
280 D0_values[current] = -1;
282 D0_colind[current] = col;
283 D0_values[current] = 1;
285 D0_rowptr[current / 2] = current;
290 LO num_coarse_edges = current / 2;
291 D0_rowptr.resize(num_coarse_edges+1);
292 D0_colind.resize(current);
293 D0_values.resize(current);
297 RCP<const Map> ownedCoarseEdgeMap = Xpetra::MapFactory<LO,GO,NO>::Build(EdgeMatrix->getRowMap()->lib(), GO_INVALID, num_coarse_edges,EdgeMatrix->getRowMap()->getIndexBase(),EdgeMatrix->getRowMap()->getComm());
301 RCP<const Map> ownedCoarseNodeMap = Pn->getDomainMap();
302 RCP<const Map> ownedPlusSharedCoarseNodeMap = D0_Pn->getCrsGraph()->getColMap();
305 RCP<CrsMatrix> D0_coarse;
310 D0_coarse = CrsMatrixFactory::Build(ownedCoarseEdgeMap,ownedPlusSharedCoarseNodeMap,0);
311 TEUCHOS_TEST_FOR_EXCEPTION(D0_coarse.is_null(),
Exceptions::RuntimeError,
"MueLu::ReitzingerPFactory: CrsMatrixFatory failed.");
317 D0_coarse->allocateAllValues(current, ia, ja, val);
318 std::copy(D0_rowptr.begin(),D0_rowptr.end(),ia.begin());
319 std::copy(D0_colind.begin(),D0_colind.end(),ja.begin());
320 std::copy(D0_values.begin(),D0_values.end(),val.begin());
321 D0_coarse->setAllValues(ia, ja, val);
326 printf(
"[%d] D0: ia.size() = %d ja.size() = %d\n",MyPID,(
int)ia.size(),(int)ja.size());
327 printf(
"[%d] D0: ia :",MyPID);
328 for(
int i=0; i<(int)ia.size(); i++)
329 printf(
"%d ",(
int)ia[i]);
330 printf(
"\n[%d] D0: global ja :",MyPID);
331 for(
int i=0; i<(int)ja.size(); i++)
332 printf(
"%d ",(
int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i]));
333 printf(
"\n[%d] D0: local ja :",MyPID);
334 for(
int i=0; i<(int)ja.size(); i++)
335 printf(
"%d ",(
int)ja[i]);
338 sprintf(fname,
"D0_global_ja_%d_%d.dat",MyPID,fineLevel.GetLevelID());
339 FILE * f = fopen(fname,
"w");
340 for(
int i=0; i<(int)ja.size(); i++)
341 fprintf(f,
"%d ",(
int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i]));
344 sprintf(fname,
"D0_local_ja_%d_%d.dat",MyPID,fineLevel.GetLevelID());
345 f = fopen(fname,
"w");
346 for(
int i=0; i<(int)ja.size(); i++)
347 fprintf(f,
"%d ",(
int)ja[i]);
354 D0_coarse->expertStaticFillComplete(ownedCoarseNodeMap,ownedCoarseEdgeMap);
356 RCP<Matrix> D0_coarse_m = rcp(
new CrsMatrixWrap(D0_coarse));
357 RCP<Teuchos::FancyOStream> fout = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
371 RCP<Matrix> Pn_D0cT = XMM::Multiply(*Pn,
false,*D0_coarse_m,
true,dummy,out0,
true,
true,
"Pn*D0c'",mm_params);
372 Pe = XMM::Multiply(*D0,
false,*Pn_D0cT,
false,dummy,out0,
true,
true,
"D0*(Pn*D0c')",mm_params);
382 SC one = Teuchos::ScalarTraits<SC>::one();
383 MT two = 2*Teuchos::ScalarTraits<MT>::one();
384 SC zero = Teuchos::ScalarTraits<SC>::zero();
387 for(LO i=0; i<(LO) Ne; i++) {
388 ArrayView<const LO> columns;
389 ArrayView<const SC> values;
390 Pe->getLocalRowView(i,columns,values);
392 ArrayView<SC> values_nc = Teuchos::av_const_cast<SC>(values);
393 for (LO j=0; j<(LO)values.size(); j++)
394 if ((values_nc[j] == one || values_nc[j] == neg_one))
399 Pe->fillComplete(Pe->getDomainMap(),Pe->getRangeMap());
403 CheckCommutingProperty(*Pe,*D0_coarse_m,*D0,*Pn);
406 Set(coarseLevel,
"P",Pe);
407 Set(coarseLevel,
"Ptent",Pe);
409 Set(coarseLevel,
"D0",D0_coarse_m);
416 int numProcs = Pe->getRowMap()->getComm()->getSize();
419 sprintf(fname,
"Pe_%d_%d.mat",numProcs,fineLevel.GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*Pe);
420 sprintf(fname,
"Pn_%d_%d.mat",numProcs,fineLevel.GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*Pn);
421 sprintf(fname,
"D0c_%d_%d.mat",numProcs,fineLevel.GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*D0_coarse_m);
422 sprintf(fname,
"D0f_%d_%d.mat",numProcs,fineLevel.GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*D0);
428 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
432 using XMM = Xpetra::MatrixMatrix<SC,LO,GO,NO>;
433 using MT =
typename Teuchos::ScalarTraits<SC>::magnitudeType;
434 SC one = Teuchos::ScalarTraits<SC>::one();
435 SC zero = Teuchos::ScalarTraits<SC>::zero();
438 Teuchos::FancyOStream &out0=GetBlackHole();
439 RCP<Matrix> left = XMM::Multiply(Pe,
false,D0_c,
false,dummy,out0);
440 RCP<Matrix> right = XMM::Multiply(D0_f,
false,Pn,
false,dummy,out0);
443 RCP<CrsMatrix> sum_c = CrsMatrixFactory::Build(left->getRowMap(),left->getNodeMaxNumRowEntries()+right->getNodeMaxNumRowEntries());
444 RCP<Matrix> summation = rcp(
new CrsMatrixWrap(sum_c));
445 XMM::TwoMatrixAdd(*left,
false, one, *summation, zero);
446 XMM::TwoMatrixAdd(*right,
false, -one, *summation, one);
448 MT norm = summation->getFrobeniusNorm();
449 GetOStream(
Statistics0) <<
"CheckCommutingProperty: ||Pe D0_c - D0_f Pn || = "<<norm<<std::endl;
461 #define MUELU_REITZINGERPFACTORY_SHORT 462 #endif // MUELU_REITZINGERPFACTORY_DEF_HPP void CheckCommutingProperty(const Matrix &Pe, const Matrix &D0_c, const Matrix &D0_f, const Matrix &Pn) const
Utility method.
Timer to be used in factories. Similar to Monitor but with additional timers.
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
Namespace for MueLu classes and methods.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
static const NoFactory * get()
Print statistics that do not involve significant additional computation.
Class that holds all level-specific information.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
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)
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.
Exception throws to report errors in the internal logical of the program.