46 #ifndef MUELU_COUPLEDAGGREGATIONCOMMHELPER_DEF_HPP 47 #define MUELU_COUPLEDAGGREGATIONCOMMHELPER_DEF_HPP 56 #include "MueLu_Utilities.hpp" 60 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
62 import_ = ImportFactory::Build(uniqueMap, nonUniqueMap);
63 tempVec_ = VectorFactory::Build(uniqueMap,
false);
66 myPID_ = uniqueMap->getComm()->getRank();
69 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 const RCP<const Map> weightMap = weight_.getMap();
72 const size_t nodeNumElements = weightMap->getNodeNumElements();
73 const RCP<const Teuchos::Comm<int> > & comm = weightMap->getComm();
74 int MyPid = comm->getRank();
78 if (comm->getSize() == 1) {
79 ArrayRCP<SC> serialWeight = weight_.getDataNonConst(0);
80 ArrayRCP<LO> serialProcWinner = procWinner_.getDataNonConst(0);
81 for (
size_t i=0; i < nodeNumElements; ++i) {
82 if (serialWeight[i] > 0) {
84 serialProcWinner[i] = MyPid;
91 #ifdef COMPARE_IN_OUT_VECTORS 92 RCP<Vector> in_weight_ = VectorFactory::Build(weight_.getMap());
94 ArrayRCP<SC> in_weight = in_weight_->getDataNonConst(0);
95 ArrayRCP<SC> weight = weight_.getDataNonConst(0);
96 for (
size_t i=0; i < nodeNumElements; ++i) in_weight[i] = weight[i];
98 RCP<LOVector> in_procWinner_ = LOVectorFactory::Build(procWinner_.getMap());
100 ArrayRCP<LO> in_procWinner = in_procWinner_->getDataNonConst(0);
101 ArrayRCP<LO> procWinner = procWinner_.getDataNonConst(0);
102 for (
size_t i=0; i < nodeNumElements; ++i) in_procWinner[i] = procWinner[i];
104 RCP<LOVector> in_companion;
106 if (companion != NULL) {
107 in_companion = LOVectorFactory::Build(companion->getMap());
108 ArrayRCP<LO> in_comp = in_companion->getDataNonConst(0);
109 ArrayRCP<LO> comp = companion->getDataNonConst(0);
110 for (
size_t i=0; i < nodeNumElements; ++i) in_comp[i] = comp[i];
115 typedef Teuchos::ScalarTraits<SC> ST;
118 if (perturbWt_ == Teuchos::null || !perturbWt_->getMap()->isSameAs(*weightMap)) {
119 perturbWt_ = VectorFactory::Build(weightMap,
false);
123 ST::seedrandom( Teuchos::as<unsigned int>(MyPid*47) );
124 for (
int i = 0; i < 10; ++i) ST::random();
129 ArrayRCP<SC> lperturbWt = perturbWt_->getDataNonConst(0);
130 for (
size_t i=0; i < nodeNumElements; ++i)
131 lperturbWt[i] = 1e-7*fabs(ST::random());
132 #ifdef COMPARE_IN_OUT_VECTORS 133 ArrayRCP<SC> locperturbWt = perturbWt_->getDataNonConst(0);
134 for (
size_t i=0; i < nodeNumElements; ++i)
135 printf(
"perturbWt[%d] = %15.10e\n",i,locperturbWt[i]);
139 ArrayRCP<SC> weight = weight_.getDataNonConst(0);
140 ArrayRCP<SC> perturbWt = perturbWt_->getDataNonConst(0);
144 SC largestGlobalWeight = weight_.normInf();
145 for (
size_t i=0; i < nodeNumElements; ++i) {
146 if (weight[i] != 0.) {
147 weight[i] += largestGlobalWeight*perturbWt[i];
159 if (postComm_ == Teuchos::null || !postComm_->getMap()->isSameAs(*weightMap) )
160 postComm_ = VectorFactory::Build(weightMap);
180 if (candidateWinners_ == Teuchos::null || !candidateWinners_->getMap()->isSameAs(*weightMap) )
181 candidateWinners_ = VectorFactory::Build(weightMap,
false);
184 ArrayRCP<SC> weight = weight_.getDataNonConst(0);
187 ArrayRCP<SC> candidateWinners = candidateWinners_->getDataNonConst(0);
188 ArrayRCP<SC> postComm = postComm_->getDataNonConst(0);
189 for (
size_t i=0; i < nodeNumElements; ++i) {
190 if (postComm[i] == weight[i]) candidateWinners[i] = (
SC) MyPid+1;
191 else candidateWinners[i] = 0;
192 weight[i]=postComm[i];
195 NonUnique2NonUnique(*candidateWinners_, *postComm_,
Xpetra::ABSMAX);
202 int numMyWinners = 0;
203 ArrayRCP<LO> procWinner = procWinner_.getDataNonConst(0);
205 ArrayRCP<SC> postComm = postComm_->getDataNonConst(0);
206 for (
size_t i=0; i < nodeNumElements; ++i) {
207 if ( weight[i] != 0.) procWinner[i] = ((int) (postComm[i])) - 1;
210 if (procWinner[i] == MyPid) ++numMyWinners;
214 weight = Teuchos::null;
216 if (companion != NULL) {
227 ArrayView<const GO> myGids = weightMap->getNodeElementList();
229 if (numMyWinners != numMyWinners_ || winnerMap_ == Teuchos::null) {
231 myWinners_ = ArrayRCP<GO>(numMyWinners);
234 numMyWinners_ = numMyWinners;
238 procWinner = Teuchos::null;
239 std::cout << MyPid <<
": nodeNumElements=" << nodeNumElements << std::endl;
240 std::cout << MyPid <<
": procWinner=" << procWinner_ << std::endl;
241 procWinner = procWinner_.getDataNonConst(0);
247 for (
size_t i = 0; i < nodeNumElements; ++i) {
248 if (procWinner[i] == MyPid) {
249 myWinners_[numMyWinners++] = myGids[i];
255 bool entryMismatch=
false;
257 for (
size_t i = 0; i < nodeNumElements; ++i) {
258 if (procWinner[i] == MyPid) {
259 if (myWinners_[numMyWinners++] != myGids[i]) {
266 if (entryMismatch ==
true) {
270 for (
size_t i = 0; i < nodeNumElements; ++i) {
271 if (procWinner[i] == MyPid) {
272 myWinners_[numMyWinners++] = myGids[i];
278 procWinner = Teuchos::null;
281 std::cout << MyPid <<
": numMyWinners=" << numMyWinners << std::endl;
282 std::cout << MyPid <<
": myWinners_" << myWinners_ << std::endl;
283 for(
int i=0;i<numMyWinners; i++)
284 std::cout << MyPid <<
": myWinners_[locId=" << i <<
"] = " << myWinners_[i] << std::endl;
290 int irealloc,orealloc;
294 if (orealloc == 1)
realloc=
true;
301 winnerMap_ = MapFactory::Build(weightMap->lib(), GSTI, myWinners_(), 0, weightMap->getComm());
307 RCP<LOVector> justWinners = LOVectorFactory::Build(winnerMap_);
310 RCP<Teuchos::FancyOStream> out = rcp(
new Teuchos::FancyOStream(rcp(&std::cout,
false)));
311 std::cout << MyPid <<
": justWinners(Vector in)=" << *justWinners << std::endl;
312 justWinners->describe(*out, Teuchos::VERB_EXTREME);
315 if ( winnerImport_ == Teuchos::null
316 || !winnerImport_->getSourceMap()->isSameAs(*weightMap)
317 || !winnerImport_->getTargetMap()->isSameAs(*winnerMap_) )
318 winnerImport_ = ImportFactory::Build(weightMap, winnerMap_);
319 RCP<const Import> winnerImport = winnerImport_;
324 catch(std::exception& e)
326 std::cout << MyPid <<
": ERR2: An exception occurred." << std::endl;
335 RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
336 fos->setOutputToRootOnly(-1);
337 if (!weightMap->getComm()->getRank())
338 std::cout <<
"------ winnerMap_ ------" << std::endl;
339 winnerMap_->describe(*fos,Teuchos::VERB_EXTREME);
340 if (!weightMap->getComm()->getRank())
341 std::cout <<
"------ weightMap ------" << std::endl;
342 weightMap->getComm()->barrier();
343 weightMap->describe(*fos,Teuchos::VERB_EXTREME);
353 if ( pushWinners_ == Teuchos::null
354 || !pushWinners_->getSourceMap()->isSameAs(*winnerMap_)
355 || !pushWinners_->getTargetMap()->isSameAs(*weightMap) )
356 pushWinners_ = ImportFactory::Build(winnerMap_,weightMap);
357 RCP<Import> pushWinners = pushWinners_;
371 catch(std::exception& e)
382 std::cout << MyPid <<
": numMyWinners=" << numMyWinners << std::endl;
384 std::cout << MyPid <<
": justWinners(Vector in)=" << std::endl;
385 justWinners->describe(*out, Teuchos::VERB_EXTREME);
387 std::cout << MyPid <<
": companion(Vector out)=" << std::endl;
388 companion->describe(*out, Teuchos::VERB_EXTREME);
391 std::cout << MyPid <<
": winnerMap_=" << *winnerMap_ << std::endl;
392 std::cout << MyPid <<
": weight_.Map=" << *weightMap << std::endl;
398 #ifdef COMPARE_IN_OUT_VECTORS 400 std::cout <<
"==============" << std::endl;
401 std::cout <<
"call #" << numCalls <<
" (1-based)" << std::endl;
402 std::cout <<
"==============" << std::endl;
409 std::string sameOrDiff;
411 ArrayRCP<SC> in_weight = in_weight_->getDataNonConst(0);
412 ArrayRCP<SC> weight = weight_.getDataNonConst(0);
413 if (MyPid == 0) std::cout <<
"==============\nweight\n==============\n" << std::endl;
414 for (
size_t i=0; i < weight_.getLocalLength(); ++i) {
415 if (in_weight[i] - weight[i] != 0) sameOrDiff =
" <<<<";
416 else sameOrDiff =
" ";
417 std::cout << std::setw(3) << i<<
": " << in_weight[i] <<
" " << weight[i] << sameOrDiff << in_weight[i] - weight[i] << std::endl;
430 ArrayRCP<LO> in_procWinner = in_procWinner_->getDataNonConst(0);
431 ArrayRCP<LO> procWinner = procWinner_.getDataNonConst(0);
432 if (MyPid == 0) std::cout <<
"==============\nprocWinner\n==============\n" << std::endl;
433 for (
size_t i=0; i < procWinner_.getLocalLength(); ++i) {
434 if (in_procWinner[i] != procWinner[i]) sameOrDiff =
" <<<<";
435 else sameOrDiff =
" ";
436 std::cout << std::setw(3) << i<<
": " << in_procWinner[i] <<
" " << procWinner[i] << sameOrDiff << std::endl;
449 if (companion != NULL) {
450 ArrayRCP<LO> in_comp = in_companion->getDataNonConst(0);
451 ArrayRCP<LO> comp = companion->getDataNonConst(0);
452 if (MyPid == 0) std::cout <<
"==============\ncompanion\n==============\n" << std::endl;
453 for (
size_t i=0; i < companion->getLocalLength(); ++i) {
454 if (in_comp[i] != comp[i]) sameOrDiff =
" <<<<";
455 else sameOrDiff =
" ";
456 std::cout << std::setw(3) << i<<
": " << in_comp[i] <<
" " << comp[i] << sameOrDiff << std::endl;
471 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
473 tempVec_->putScalar(0.);
477 tempVec_->doExport(source, *import_, what);
480 catch(std::exception& e)
482 int MyPid = tempVec_->getMap()->getComm()->getRank();
483 std::cout << MyPid <<
": ERR1: An exception occurred." << std::endl;
490 #endif // MUELU_COUPLEDAGGREGATIONCOMMHELPER_DEF_HPP
#define MueLu_maxAll(rcpComm, in, out)
Namespace for MueLu classes and methods.
void NonUnique2NonUnique(const Vector &source, Vector &dest, const Xpetra::CombineMode what) const
Redistribute data in source to dest where both source and dest might have multiple copies of the same...
CoupledAggregationCommHelper(const RCP< const Map > &uniqueMap, const RCP< const Map > &nonUniqueMap)
Constructor.
void realloc(DynRankView< T, P... > &v, const size_t n0=~size_t(0), const size_t n1=~size_t(0), const size_t n2=~size_t(0), const size_t n3=~size_t(0), const size_t n4=~size_t(0), const size_t n5=~size_t(0), const size_t n6=~size_t(0), const size_t n7=~size_t(0))
void ArbitrateAndCommunicate(Vector &weights, Aggregates &aggregates, const bool perturb) const
This method assigns unknowns to aggregates.