46 #ifndef MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP 47 #define MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP 49 #include <Xpetra_Map.hpp> 50 #include <Xpetra_MultiVectorFactory.hpp> 51 #include <Xpetra_VectorFactory.hpp> 58 #include "MueLu_CoupledAggregationCommHelper.hpp" 64 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
66 phase3AggCreation_(.5),
67 minNodesPerAggregate_(1)
70 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 Monitor m(*
this,
"AggregateLeftovers");
75 int exp_nRows = aggregates.
GetMap()->getNodeNumElements();
76 int myPid = graph.
GetComm()->getRank();
79 int minNodesPerAggregate = GetMinNodesPerAggregate();
81 const RCP<const Map> nonUniqueMap = aggregates.
GetMap();
87 RCP<Xpetra::Vector<double,LO,GO,NO> > distWeights = Xpetra::VectorFactory<double,LO,GO,NO>::Build(nonUniqueMap);
93 ArrayRCP<const LO> vertex2AggId = aggregates.
GetVertex2AggId()->getData(0);
94 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
95 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
97 for (
size_t i=0;i<nonUniqueMap->getNodeNumElements();i++) {
101 if (aggregates.
IsRoot(i)) weights[i] = 2.;
116 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
117 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
118 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
120 for (
my_size_t i = 0; i < nVertices; i++) {
121 if ( aggregates.
IsRoot(i) && (procWinner[i] == myPid) ) {
126 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
130 vertex2AggId[colj] = vertex2AggId[i];
144 GO total_phase_one_aggregated = 0;
146 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
148 GO phase_one_aggregated = 0;
149 for (
my_size_t i = 0; i < nVertices; i++) {
151 phase_one_aggregated++;
156 GO local_nVertices = nVertices, total_nVertices = 0;
165 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
168 factor = ((double) total_phase_one_aggregated)/((double)(total_nVertices + 1));
169 factor = pow(factor, GetPhase3AggCreation());
171 for (
my_size_t i = 0; i < nVertices; i++) {
177 int rowi_N = neighOfINode.size();
179 int nonaggd_neighbors = 0;
180 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
185 if ( (nonaggd_neighbors > minNodesPerAggregate) &&
186 (((double) nonaggd_neighbors)/((double) rowi_N) > factor))
188 vertex2AggId[i] = (nAggregates)++;
189 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
192 vertex2AggId[colj] = vertex2AggId[i];
193 if (colj < nVertices) weights[colj] = 2.;
194 else weights[colj] = 1.;
211 GO Nphase1_agg = nAggregates;
216 GetOStream(
Statistics1) <<
"Phase 1 - nodes aggregated = " << total_phase_one_aggregated << std::endl;
217 GetOStream(
Statistics1) <<
"Phase 1 - total aggregates = " << total_aggs << std::endl;
219 GO i = nAggregates - Nphase1_agg;
221 GetOStream(
Statistics1) <<
"Phase 3 - additional aggregates = " << i << std::endl;
231 RCP<Xpetra::Vector<double,LO,GO,NO> > temp_ = Xpetra::VectorFactory<double,LO,GO,NO> ::Build(nonUniqueMap,
false);
232 temp_->putScalar(1.);
234 RCP<Xpetra::Vector<double,LO,GO,NO> > tempOutput_ = Xpetra::VectorFactory<double,LO,GO,NO> ::Build(nonUniqueMap);
238 std::vector<bool> gidNotShared(exp_nRows);
240 ArrayRCP<const double> tempOutput = tempOutput_->getData(0);
241 for (
int i = 0; i < exp_nRows; i++) {
242 if (tempOutput[i] > 1.)
243 gidNotShared[i] =
false;
245 gidNotShared[i] =
true;
250 double nAggregatesTarget;
251 nAggregatesTarget = ((double) uniqueMap->getGlobalNumElements())* (((
double) uniqueMap->getGlobalNumElements())/ ((
double) graph.
GetGlobalNumEdges()));
253 GO nAggregatesLocal=nAggregates, nAggregatesGlobal;
MueLu_sumAll(graph.
GetComm(), nAggregatesLocal, nAggregatesGlobal);
262 #define MUELU_PHASE4BUCKETS 6 263 if ((nAggregatesGlobal < graph.
GetComm()->getSize()) &&
264 (2.5*nAggregatesGlobal < nAggregatesTarget) &&
265 (minNAggs ==0) && (maxNAggs <= 1)) {
269 typedef Teuchos::ScalarTraits<double> scalarTrait;
270 scalarTrait::seedrandom(static_cast<unsigned int>(myPid*2 + (
int) (11*scalarTrait::random())));
271 int k = (int)ceil( (10.*myPid)/graph.
GetComm()->getSize());
272 for (
int i = 0; i < k+7; i++) scalarTrait::random();
273 temp_->setSeed(static_cast<unsigned int>(scalarTrait::random()));
278 ArrayRCP<double> temp = temp_->getDataNonConst(0);
286 ArrayRCP<LO> candidates = Teuchos::arcp<LO>(nVertices+1);
288 double priorThreshold = 0.;
292 ArrayRCP<const LO> vertex2AggId = aggregates.
GetVertex2AggId()->getData(0);
293 ArrayView<const LO> vertex2AggIdView = vertex2AggId();
294 RootCandidates(nVertices, vertex2AggIdView, graph, candidates, nCandidates, nCandidatesGlobal);
298 double nTargetNewGuys = nAggregatesTarget - nAggregatesGlobal;
299 double threshold = priorThreshold + (1. - priorThreshold)*nTargetNewGuys/(nCandidatesGlobal + .001);
302 priorThreshold = threshold;
305 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
306 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
308 for (
int k = 0; k < nCandidates; k++ ) {
309 int i = candidates[k];
316 if (neighOfINode.size() > minNodesPerAggregate) {
318 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
324 vertex2AggId[Adjacent] = nAggregates;
325 weights[Adjacent] = 1.;
328 if (count >= minNodesPerAggregate) {
329 vertex2AggId[i] = nAggregates++;
334 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
336 if (vertex2AggId[Adjacent] == nAggregates){
338 weights[Adjacent] = 0.;
350 nAggregatesLocal=nAggregates;
357 RemoveSmallAggs(aggregates, minNodesPerAggregate, distWeights, myWidget);
372 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
373 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
374 for(LO k = 0; k < vertex2AggId.size(); ++k )
375 if(vertex2AggId[k]>observedNAgg)
376 observedNAgg=vertex2AggId[k];
380 ArrayRCP<int> Mark = Teuchos::arcp<int>(exp_nRows+1);
381 ArrayRCP<int> agg_incremented = Teuchos::arcp<int>(observedNAgg);
382 ArrayRCP<int> SumOfMarks = Teuchos::arcp<int>(observedNAgg);
385 for (
int i = 0; i < agg_incremented.size(); i++) agg_incremented[i] = 0;
386 for (
int i = 0; i < SumOfMarks.size(); i++) SumOfMarks[i] = 0;
390 std::vector<int> RowPtr(exp_nRows+1-nVertices);
392 ArrayRCP<const LO> vertex2AggIdCst = aggregates.
GetVertex2AggId()->getData(0);
394 for (
int i = nVertices; i < exp_nRows; i++) RowPtr[i-nVertices] = 0;
395 for (
int i = 0; i < nVertices; i++) {
400 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
403 RowPtr[j-nVertices]++;
411 for (
int i = nVertices; i < exp_nRows; i++) {
412 iTemp = RowPtr[i-nVertices];
413 RowPtr[i-nVertices] = iSum;
416 RowPtr[exp_nRows-nVertices] = iSum;
417 std::vector<LO> cols(iSum+1);
420 for (
int i = 0; i < nVertices; i++) {
425 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
428 cols[RowPtr[j-nVertices]++] = i;
434 for (
int i = exp_nRows; i > nVertices; i--)
435 RowPtr[i-nVertices] = RowPtr[i-1-nVertices];
439 vertex2AggIdCst = Teuchos::null;
443 int thresholds[10] = {300,200,100,50,25,13,7,4,2,0};
450 for (
int kk = 0; kk < 10; kk += 2) {
451 bestScoreCutoff = thresholds[kk];
453 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
454 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
455 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
457 for (
int i = 0; i < exp_nRows; i++) {
462 ArrayView<const LO> neighOfINode;
470 LO *rowi_col = NULL, rowi_N;
471 rowi_col = &(cols[RowPtr[i-nVertices]]);
472 rowi_N = RowPtr[i+1-nVertices] - RowPtr[i-nVertices];
474 neighOfINode = ArrayView<const LO>(rowi_col, rowi_N);
476 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
478 int AdjacentAgg = vertex2AggId[Adjacent];
483 ((procWinner[Adjacent] == myPid) ||
485 SumOfMarks[AdjacentAgg] += Mark[Adjacent];
491 bool cannotLoseAllFriends=
false;
493 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
495 int AdjacentAgg = vertex2AggId[Adjacent];
499 (SumOfMarks[AdjacentAgg] != 0) &&
500 ((procWinner[Adjacent] == myPid) ||
507 double penalty = (double) (
INCR_SCALING*agg_incremented[AdjacentAgg]);
510 int score = SumOfMarks[AdjacentAgg]- ((int) floor(penalty));
512 if (score > best_score) {
513 best_agg = AdjacentAgg;
515 BestMark = Mark[Adjacent];
516 cannotLoseAllFriends =
false;
523 if ((weights[Adjacent]== 0.) || (gidNotShared[Adjacent] ==
true))
524 cannotLoseAllFriends =
true;
529 else if (best_agg == AdjacentAgg) {
530 if ((weights[Adjacent]== 0.) || (gidNotShared[Adjacent] ==
true))
531 cannotLoseAllFriends =
true;
532 if (Mark[Adjacent] > BestMark) BestMark = Mark[Adjacent];
537 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
539 int AdjacentAgg = vertex2AggId[Adjacent];
540 if (AdjacentAgg >= 0) SumOfMarks[AdjacentAgg] = 0;
543 if ( (best_score >= bestScoreCutoff) && (cannotLoseAllFriends)) {
547 vertex2AggId[i] = best_agg;
548 weights[i] = best_score;
549 agg_incremented[best_agg]++;
550 Mark[i] = (int) ceil( ((
double) BestMark)/2.);
557 vertex2AggId = Teuchos::null;
558 procWinner = Teuchos::null;
559 weights = Teuchos::null;
579 int Nleftover = 0, Nsingle = 0;
582 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
583 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
584 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
587 for (
my_size_t i = 0; i < nVertices; i++) {
597 vertex2AggId[i] = nAggregates;
601 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
605 vertex2AggId[j] = nAggregates;
610 if ( count >= minNodesPerAggregate) {
621 if (nAggregates > 0) {
622 for (
my_size_t i = 0; i < nVertices; i++) {
623 if ((vertex2AggId[i] == nAggregates) && (procWinner[i] == myPid)) {
635 if (nAggregates > 0) {
636 for (
my_size_t i = 0; i < nVertices; i++) {
643 if (vertex2AggId[i] == nAggregates) {
644 vertex2AggId[i] = nAggregates-1;
667 GetOStream(
Statistics1) <<
"Phase 6 - leftovers = " << total_Nleftover <<
" and singletons = " << total_Nsingle << std::endl;
674 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
676 ArrayView<const LO> & vertex2AggId,
GraphBase const &graph,
681 for (
my_size_t i = 0; i < nVertices; i++ ) {
683 bool noAggdNeighbors =
true;
688 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
691 noAggdNeighbors =
false;
693 if (noAggdNeighbors ==
true) candidates[nCandidates++] = i;
701 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
704 int myPid = aggregates.
GetMap()->getComm()->getRank();
708 ArrayRCP<LO> procWinner = aggregates.
GetProcWinner()->getDataNonConst(0);
709 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
710 LO size = procWinner.size();
716 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
722 for (LO i = 0; i < nAggregates; i++) {
723 if ( AggInfo[i] < min_size) {
726 else AggInfo[i] = NewNAggs++;
729 for (LO k = 0; k < size; k++ ) {
730 if (procWinner[k] == myPid) {
732 vertex2AggId[k] = AggInfo[vertex2AggId[k]];
739 nAggregates = NewNAggs;
747 for (LO i = 0; i < size; i++) {
758 #endif // MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP
#define MueLu_sumAll(rcpComm, in, out)
#define MueLu_maxAll(rcpComm, in, out)
Container class for aggregation information.
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
bool IsRoot(LO i) const
Returns true if node with given local node id is marked to be a root node.
virtual const RCP< const Map > GetDomainMap() const =0
Namespace for MueLu classes and methods.
virtual Xpetra::global_size_t GetGlobalNumEdges() const =0
Return number of global edges in the graph.
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...
#define MueLu_minAll(rcpComm, in, out)
void SetIsRoot(LO i, bool value=true)
Set root node information.
Helper class for providing arbitrated communication across processors.
#define MUELU_UNAGGREGATED
virtual const RCP< const Teuchos::Comm< int > > GetComm() const =0
const RCP< const Map > GetMap() const
returns (overlapping) map of aggregate/node distribution
MueLu representation of a graph.
void AggregateLeftovers(GraphBase const &graph, Aggregates &aggregates) const
Take a partially aggregated graph and complete the aggregation.
LeftoverAggregationAlgorithm()
Constructor.
Timer to be used in non-factories.
LO GetNumAggregates() const
returns the number of aggregates of the current processor. Note: could/should be renamed to GetNumLoc...
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
#define MUELU_PHASE4BUCKETS
#define MUELU_DISTONE_VERTEX_WEIGHT
Exception throws to report errors in the internal logical of the program.
const RCP< LOVector > & GetProcWinner() const
Returns constant vector that maps local node IDs to owning processor IDs.
#define MUELU_PENALTYFACTOR
int RemoveSmallAggs(Aggregates &aggregates, int min_size, RCP< Xpetra::Vector< double, LO, GO, NO > > &distWeights, const MueLu::CoupledAggregationCommHelper< LO, GO, NO > &myWidget) const
Attempt to clean up aggregates that are too small.
void ArbitrateAndCommunicate(Vector &weights, Aggregates &aggregates, const bool perturb) const
This method assigns unknowns to aggregates.
void RootCandidates(my_size_t nVertices, ArrayView< const LO > &vertex2AggId, GraphBase const &graph, ArrayRCP< LO > &candidates, my_size_t &nCandidates, global_size_t &nCandidatesGlobal) const
Build a list of candidate root nodes.
virtual Teuchos::ArrayView< const LocalOrdinal > getNeighborVertices(LocalOrdinal v) const =0
Return the list of vertices adjacent to the vertex 'v'.
Teuchos::ArrayRCP< LO > ComputeAggregateSizes(bool forceRecompute=false) const
Compute sizes of aggregates.
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.