46 #ifndef MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP 47 #define MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP 57 #include "MueLu_CoupledAggregationCommHelper.hpp" 63 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
65 phase3AggCreation_(.5),
66 minNodesPerAggregate_(1)
69 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 Monitor m(*
this,
"AggregateLeftovers");
74 int exp_nRows = aggregates.
GetMap()->getNodeNumElements();
75 int myPid = graph.
GetComm()->getRank();
78 int minNodesPerAggregate = GetMinNodesPerAggregate();
80 const RCP<const Map> nonUniqueMap = aggregates.
GetMap();
92 ArrayRCP<const LO> vertex2AggId = aggregates.
GetVertex2AggId()->getData(0);
93 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
94 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
96 for (
size_t i=0;i<nonUniqueMap->getNodeNumElements();i++) {
100 if (aggregates.
IsRoot(i)) weights[i] = 2.;
115 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
116 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
117 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
119 for (
my_size_t i = 0; i < nVertices; i++) {
120 if ( aggregates.
IsRoot(i) && (procWinner[i] == myPid) ) {
125 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
129 vertex2AggId[colj] = vertex2AggId[i];
143 GO total_phase_one_aggregated = 0;
145 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
147 GO phase_one_aggregated = 0;
148 for (
my_size_t i = 0; i < nVertices; i++) {
150 phase_one_aggregated++;
155 GO local_nVertices = nVertices, total_nVertices = 0;
164 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
167 factor = ((double) total_phase_one_aggregated)/((double)(total_nVertices + 1));
168 factor = pow(factor, GetPhase3AggCreation());
170 for (
my_size_t i = 0; i < nVertices; i++) {
176 int rowi_N = neighOfINode.size();
178 int nonaggd_neighbors = 0;
179 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
184 if ( (nonaggd_neighbors > minNodesPerAggregate) &&
185 (((double) nonaggd_neighbors)/((double) rowi_N) > factor))
187 vertex2AggId[i] = (nAggregates)++;
188 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
191 vertex2AggId[colj] = vertex2AggId[i];
192 if (colj < nVertices) weights[colj] = 2.;
193 else weights[colj] = 1.;
210 GO Nphase1_agg = nAggregates;
215 GetOStream(
Statistics1) <<
"Phase 1 - nodes aggregated = " << total_phase_one_aggregated << std::endl;
216 GetOStream(
Statistics1) <<
"Phase 1 - total aggregates = " << total_aggs << std::endl;
218 GO i = nAggregates - Nphase1_agg;
220 GetOStream(
Statistics1) <<
"Phase 3 - additional aggregates = " << i << std::endl;
231 temp_->putScalar(1.);
237 std::vector<bool> gidNotShared(exp_nRows);
239 ArrayRCP<const double> tempOutput = tempOutput_->getData(0);
240 for (
int i = 0; i < exp_nRows; i++) {
241 if (tempOutput[i] > 1.)
242 gidNotShared[i] =
false;
244 gidNotShared[i] =
true;
249 double nAggregatesTarget;
250 nAggregatesTarget = ((double) uniqueMap->getGlobalNumElements())* (((
double) uniqueMap->getGlobalNumElements())/ ((
double) graph.
GetGlobalNumEdges()));
252 GO nAggregatesLocal=nAggregates, nAggregatesGlobal;
MueLu_sumAll(graph.
GetComm(), nAggregatesLocal, nAggregatesGlobal);
261 #define MUELU_PHASE4BUCKETS 6 262 if ((nAggregatesGlobal < graph.
GetComm()->getSize()) &&
263 (2.5*nAggregatesGlobal < nAggregatesTarget) &&
264 (minNAggs ==0) && (maxNAggs <= 1)) {
268 typedef Teuchos::ScalarTraits<double> scalarTrait;
269 scalarTrait::seedrandom(static_cast<unsigned int>(myPid*2 + (
int) (11*scalarTrait::random())));
270 int k = (int)ceil( (10.*myPid)/graph.
GetComm()->getSize());
271 for (
int i = 0; i < k+7; i++) scalarTrait::random();
272 temp_->setSeed(static_cast<unsigned int>(scalarTrait::random()));
277 ArrayRCP<double> temp = temp_->getDataNonConst(0);
285 ArrayRCP<LO> candidates = Teuchos::arcp<LO>(nVertices+1);
287 double priorThreshold = 0.;
291 ArrayRCP<const LO> vertex2AggId = aggregates.
GetVertex2AggId()->getData(0);
292 ArrayView<const LO> vertex2AggIdView = vertex2AggId();
293 RootCandidates(nVertices, vertex2AggIdView, graph, candidates, nCandidates, nCandidatesGlobal);
297 double nTargetNewGuys = nAggregatesTarget - nAggregatesGlobal;
298 double threshold = priorThreshold + (1. - priorThreshold)*nTargetNewGuys/(nCandidatesGlobal + .001);
301 priorThreshold = threshold;
304 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
305 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
307 for (
int k = 0; k < nCandidates; k++ ) {
308 int i = candidates[k];
315 if (neighOfINode.size() > minNodesPerAggregate) {
317 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
323 vertex2AggId[Adjacent] = nAggregates;
324 weights[Adjacent] = 1.;
327 if (count >= minNodesPerAggregate) {
328 vertex2AggId[i] = nAggregates++;
333 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
335 if (vertex2AggId[Adjacent] == nAggregates){
337 weights[Adjacent] = 0.;
349 nAggregatesLocal=nAggregates;
356 RemoveSmallAggs(aggregates, minNodesPerAggregate, distWeights, myWidget);
371 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
372 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
373 for(
LO k = 0; k < vertex2AggId.size(); ++k )
374 if(vertex2AggId[k]>observedNAgg)
375 observedNAgg=vertex2AggId[k];
379 ArrayRCP<int> Mark = Teuchos::arcp<int>(exp_nRows+1);
380 ArrayRCP<int> agg_incremented = Teuchos::arcp<int>(observedNAgg);
381 ArrayRCP<int> SumOfMarks = Teuchos::arcp<int>(observedNAgg);
384 for (
int i = 0; i < agg_incremented.size(); i++) agg_incremented[i] = 0;
385 for (
int i = 0; i < SumOfMarks.size(); i++) SumOfMarks[i] = 0;
389 std::vector<int> RowPtr(exp_nRows+1-nVertices);
391 ArrayRCP<const LO> vertex2AggIdCst = aggregates.
GetVertex2AggId()->getData(0);
393 for (
int i = nVertices; i < exp_nRows; i++) RowPtr[i-nVertices] = 0;
394 for (
int i = 0; i < nVertices; i++) {
399 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
402 RowPtr[j-nVertices]++;
410 for (
int i = nVertices; i < exp_nRows; i++) {
411 iTemp = RowPtr[i-nVertices];
412 RowPtr[i-nVertices] = iSum;
415 RowPtr[exp_nRows-nVertices] = iSum;
416 std::vector<LO> cols(iSum+1);
419 for (
int i = 0; i < nVertices; i++) {
424 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
427 cols[RowPtr[j-nVertices]++] = i;
433 for (
int i = exp_nRows; i > nVertices; i--)
434 RowPtr[i-nVertices] = RowPtr[i-1-nVertices];
438 vertex2AggIdCst = Teuchos::null;
442 int thresholds[10] = {300,200,100,50,25,13,7,4,2,0};
449 for (
int kk = 0; kk < 10; kk += 2) {
450 bestScoreCutoff = thresholds[kk];
452 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
453 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
454 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
456 for (
int i = 0; i < exp_nRows; i++) {
461 ArrayView<const LO> neighOfINode;
469 LO *rowi_col = NULL, rowi_N;
470 rowi_col = &(cols[RowPtr[i-nVertices]]);
471 rowi_N = RowPtr[i+1-nVertices] - RowPtr[i-nVertices];
473 neighOfINode = ArrayView<const LO>(rowi_col, rowi_N);
475 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
477 int AdjacentAgg = vertex2AggId[Adjacent];
482 ((procWinner[Adjacent] == myPid) ||
484 SumOfMarks[AdjacentAgg] += Mark[Adjacent];
490 bool cannotLoseAllFriends=
false;
492 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
494 int AdjacentAgg = vertex2AggId[Adjacent];
498 (SumOfMarks[AdjacentAgg] != 0) &&
499 ((procWinner[Adjacent] == myPid) ||
506 double penalty = (double) (
INCR_SCALING*agg_incremented[AdjacentAgg]);
509 int score = SumOfMarks[AdjacentAgg]- ((int) floor(penalty));
511 if (score > best_score) {
512 best_agg = AdjacentAgg;
514 BestMark = Mark[Adjacent];
515 cannotLoseAllFriends =
false;
522 if ((weights[Adjacent]== 0.) || (gidNotShared[Adjacent] ==
true))
523 cannotLoseAllFriends =
true;
528 else if (best_agg == AdjacentAgg) {
529 if ((weights[Adjacent]== 0.) || (gidNotShared[Adjacent] ==
true))
530 cannotLoseAllFriends =
true;
531 if (Mark[Adjacent] > BestMark) BestMark = Mark[Adjacent];
536 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
538 int AdjacentAgg = vertex2AggId[Adjacent];
539 if (AdjacentAgg >= 0) SumOfMarks[AdjacentAgg] = 0;
542 if ( (best_score >= bestScoreCutoff) && (cannotLoseAllFriends)) {
546 vertex2AggId[i] = best_agg;
547 weights[i] = best_score;
548 agg_incremented[best_agg]++;
549 Mark[i] = (int) ceil( ((
double) BestMark)/2.);
556 vertex2AggId = Teuchos::null;
557 procWinner = Teuchos::null;
558 weights = Teuchos::null;
578 int Nleftover = 0, Nsingle = 0;
581 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
582 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
583 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
586 for (
my_size_t i = 0; i < nVertices; i++) {
596 vertex2AggId[i] = nAggregates;
600 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
604 vertex2AggId[j] = nAggregates;
609 if ( count >= minNodesPerAggregate) {
620 if (nAggregates > 0) {
621 for (
my_size_t i = 0; i < nVertices; i++) {
622 if ((vertex2AggId[i] == nAggregates) && (procWinner[i] == myPid)) {
634 if (nAggregates > 0) {
635 for (
my_size_t i = 0; i < nVertices; i++) {
642 if (vertex2AggId[i] == nAggregates) {
643 vertex2AggId[i] = nAggregates-1;
666 GetOStream(
Statistics1) <<
"Phase 6 - leftovers = " << total_Nleftover <<
" and singletons = " << total_Nsingle << std::endl;
673 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
675 ArrayView<const LO> & vertex2AggId,
GraphBase const &graph,
680 for (
my_size_t i = 0; i < nVertices; i++ ) {
682 bool noAggdNeighbors =
true;
687 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
690 noAggdNeighbors =
false;
692 if (noAggdNeighbors ==
true) candidates[nCandidates++] = i;
700 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
703 int myPid = aggregates.
GetMap()->getComm()->getRank();
707 ArrayRCP<LO> procWinner = aggregates.
GetProcWinner()->getDataNonConst(0);
708 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
709 LO size = procWinner.size();
715 ArrayRCP<double> weights = distWeights->getDataNonConst(0);
721 for (
LO i = 0; i < nAggregates; i++) {
722 if ( AggInfo[i] < min_size) {
725 else AggInfo[i] = NewNAggs++;
728 for (
LO k = 0; k < size; k++ ) {
729 if (procWinner[k] == myPid) {
731 vertex2AggId[k] = AggInfo[vertex2AggId[k]];
738 nAggregates = NewNAggs;
746 for (
LO i = 0; i < size; i++) {
757 #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.
const RCP< LOVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
#define MUELU_UNAGGREGATED
virtual const RCP< const Teuchos::Comm< int > > GetComm() const =0
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
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...
#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'.
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.
Teuchos::ArrayRCP< LO > ComputeAggregateSizes(bool forceRecompute=true, bool cacheSizes=false) const
Compute sizes of aggregates.