46 #ifndef MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP
47 #define MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP
49 #include <Xpetra_Map.hpp>
50 #include <Xpetra_MultiVectorFactory.hpp>
58 #include "MueLu_CoupledAggregationCommHelper.hpp"
62 #ifdef HAVE_MUELU_TPETRA
63 #include <Tpetra_Details_DefaultTypes.hpp>
69 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 phase3AggCreation_(.5),
72 minNodesPerAggregate_(1)
75 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
77 Monitor m(*
this,
"AggregateLeftovers");
80 int exp_nRows = aggregates.
GetMap()->getLocalNumElements();
81 int myPid = graph.
GetComm()->getRank();
84 int minNodesPerAggregate = GetMinNodesPerAggregate();
102 for (
size_t i=0;i<nonUniqueMap->getLocalNumElements();i++) {
106 if (aggregates.
IsRoot(i)) weights[i] = 2.;
125 for (
my_size_t i = 0; i < nVertices; i++) {
126 if ( aggregates.
IsRoot(i) && (procWinner[i] == myPid) ) {
135 vertex2AggId[colj] = vertex2AggId[i];
149 GO total_phase_one_aggregated = 0;
153 GO phase_one_aggregated = 0;
154 for (
my_size_t i = 0; i < nVertices; i++) {
156 phase_one_aggregated++;
161 GO local_nVertices = nVertices, total_nVertices = 0;
173 factor = ((double) total_phase_one_aggregated)/((double)(total_nVertices + 1));
174 factor = pow(factor, GetPhase3AggCreation());
176 for (
my_size_t i = 0; i < nVertices; i++) {
182 int rowi_N = neighOfINode.
size();
184 int nonaggd_neighbors = 0;
190 if ( (nonaggd_neighbors > minNodesPerAggregate) &&
191 (((double) nonaggd_neighbors)/((double) rowi_N) > factor))
193 vertex2AggId[i] = (nAggregates)++;
197 vertex2AggId[colj] = vertex2AggId[i];
198 if (colj < nVertices) weights[colj] = 2.;
199 else weights[colj] = 1.;
216 GO Nphase1_agg = nAggregates;
221 GetOStream(
Statistics1) <<
"Phase 1 - nodes aggregated = " << total_phase_one_aggregated << std::endl;
222 GetOStream(
Statistics1) <<
"Phase 1 - total aggregates = " << total_aggs << std::endl;
224 GO i = nAggregates - Nphase1_agg;
226 GetOStream(
Statistics1) <<
"Phase 3 - additional aggregates = " << i << std::endl;
237 temp_->putScalar(1.);
243 std::vector<bool> gidNotShared(exp_nRows);
246 for (
int i = 0; i < exp_nRows; i++) {
247 if (tempOutput[i] > 1.)
248 gidNotShared[i] =
false;
250 gidNotShared[i] =
true;
255 double nAggregatesTarget;
256 nAggregatesTarget = ((double) uniqueMap->getGlobalNumElements())* (((
double) uniqueMap->getGlobalNumElements())/ ((
double) graph.
GetGlobalNumEdges()));
258 GO nAggregatesLocal=nAggregates, nAggregatesGlobal;
MueLu_sumAll(graph.
GetComm(), nAggregatesLocal, nAggregatesGlobal);
267 #define MUELU_PHASE4BUCKETS 6
268 if ((nAggregatesGlobal < graph.
GetComm()->getSize()) &&
269 (2.5*nAggregatesGlobal < nAggregatesTarget) &&
270 (minNAggs ==0) && (maxNAggs <= 1)) {
275 scalarTrait::seedrandom(static_cast<unsigned int>(myPid*2 + (
int) (11*scalarTrait::random())));
276 int k = (int)ceil( (10.*myPid)/graph.
GetComm()->getSize());
277 for (
int i = 0; i < k+7; i++) scalarTrait::random();
278 temp_->setSeed(static_cast<unsigned int>(scalarTrait::random()));
291 ArrayRCP<LO> candidates = Teuchos::arcp<LO>(nVertices+1);
293 double priorThreshold = 0.;
299 RootCandidates(nVertices, vertex2AggIdView, graph, candidates, nCandidates, nCandidatesGlobal);
303 double nTargetNewGuys = nAggregatesTarget - nAggregatesGlobal;
304 double threshold = priorThreshold + (1. - priorThreshold)*nTargetNewGuys/(nCandidatesGlobal + .001);
307 priorThreshold = threshold;
313 for (
int k = 0; k < nCandidates; k++ ) {
314 int i = candidates[k];
321 if (neighOfINode.
size() > minNodesPerAggregate) {
329 vertex2AggId[Adjacent] = nAggregates;
330 weights[Adjacent] = 1.;
333 if (count >= minNodesPerAggregate) {
334 vertex2AggId[i] = nAggregates++;
341 if (vertex2AggId[Adjacent] == nAggregates){
343 weights[Adjacent] = 0.;
355 nAggregatesLocal=nAggregates;
362 RemoveSmallAggs(aggregates, minNodesPerAggregate, distWeights, myWidget);
379 for(
LO k = 0; k < vertex2AggId.
size(); ++k )
380 if(vertex2AggId[k]>observedNAgg)
381 observedNAgg=vertex2AggId[k];
386 ArrayRCP<int> agg_incremented = Teuchos::arcp<int>(observedNAgg);
390 for (
int i = 0; i < agg_incremented.size(); i++) agg_incremented[i] = 0;
391 for (
int i = 0; i < SumOfMarks.size(); i++) SumOfMarks[i] = 0;
395 std::vector<int> RowPtr(exp_nRows+1-nVertices);
399 for (
int i = nVertices; i < exp_nRows; i++) RowPtr[i-nVertices] = 0;
400 for (
int i = 0; i < nVertices; i++) {
408 RowPtr[j-nVertices]++;
416 for (
int i = nVertices; i < exp_nRows; i++) {
417 iTemp = RowPtr[i-nVertices];
418 RowPtr[i-nVertices] = iSum;
421 RowPtr[exp_nRows-nVertices] = iSum;
422 std::vector<LO> cols(iSum+1);
425 for (
int i = 0; i < nVertices; i++) {
433 cols[RowPtr[j-nVertices]++] = i;
439 for (
int i = exp_nRows; i > nVertices; i--)
440 RowPtr[i-nVertices] = RowPtr[i-1-nVertices];
444 vertex2AggIdCst = Teuchos::null;
448 int thresholds[10] = {300,200,100,50,25,13,7,4,2,0};
455 for (
int kk = 0; kk < 10; kk += 2) {
456 bestScoreCutoff = thresholds[kk];
462 for (
int i = 0; i < exp_nRows; i++) {
475 LO *rowi_col = NULL, rowi_N;
476 rowi_col = &(cols[RowPtr[i-nVertices]]);
477 rowi_N = RowPtr[i+1-nVertices] - RowPtr[i-nVertices];
483 int AdjacentAgg = vertex2AggId[Adjacent];
488 ((procWinner[Adjacent] == myPid) ||
490 SumOfMarks[AdjacentAgg] += Mark[Adjacent];
496 bool cannotLoseAllFriends=
false;
500 int AdjacentAgg = vertex2AggId[Adjacent];
504 (SumOfMarks[AdjacentAgg] != 0) &&
505 ((procWinner[Adjacent] == myPid) ||
512 double penalty = (double) (
INCR_SCALING*agg_incremented[AdjacentAgg]);
515 int score = SumOfMarks[AdjacentAgg]- ((int) floor(penalty));
517 if (score > best_score) {
518 best_agg = AdjacentAgg;
520 BestMark = Mark[Adjacent];
521 cannotLoseAllFriends =
false;
528 if ((weights[Adjacent]== 0.) || (gidNotShared[Adjacent] ==
true))
529 cannotLoseAllFriends =
true;
534 else if (best_agg == AdjacentAgg) {
535 if ((weights[Adjacent]== 0.) || (gidNotShared[Adjacent] ==
true))
536 cannotLoseAllFriends =
true;
537 if (Mark[Adjacent] > BestMark) BestMark = Mark[Adjacent];
544 int AdjacentAgg = vertex2AggId[Adjacent];
545 if (AdjacentAgg >= 0) SumOfMarks[AdjacentAgg] = 0;
548 if ( (best_score >= bestScoreCutoff) && (cannotLoseAllFriends)) {
552 vertex2AggId[i] = best_agg;
553 weights[i] = best_score;
554 agg_incremented[best_agg]++;
555 Mark[i] = (int) ceil( ((
double) BestMark)/2.);
562 vertex2AggId = Teuchos::null;
563 procWinner = Teuchos::null;
564 weights = Teuchos::null;
584 int Nleftover = 0, Nsingle = 0;
592 for (
my_size_t i = 0; i < nVertices; i++) {
602 vertex2AggId[i] = nAggregates;
610 vertex2AggId[j] = nAggregates;
615 if ( count >= minNodesPerAggregate) {
626 if (nAggregates > 0) {
627 for (
my_size_t i = 0; i < nVertices; i++) {
628 if ((vertex2AggId[i] == nAggregates) && (procWinner[i] == myPid)) {
640 if (nAggregates > 0) {
641 for (
my_size_t i = 0; i < nVertices; i++) {
648 if (vertex2AggId[i] == nAggregates) {
649 vertex2AggId[i] = nAggregates-1;
672 GetOStream(
Statistics1) <<
"Phase 6 - leftovers = " << total_Nleftover <<
" and singletons = " << total_Nsingle << std::endl;
679 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
686 for (
my_size_t i = 0; i < nVertices; i++ ) {
688 bool noAggdNeighbors =
true;
696 noAggdNeighbors =
false;
698 if (noAggdNeighbors ==
true) candidates[nCandidates++] = i;
706 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
709 int myPid = aggregates.
GetMap()->getComm()->getRank();
715 LO size = procWinner.
size();
727 for (
LO i = 0; i < nAggregates; i++) {
728 if ( AggInfo[i] < min_size) {
731 else AggInfo[i] = NewNAggs++;
734 for (
LO k = 0; k < size; k++ ) {
735 if (procWinner[k] == myPid) {
737 vertex2AggId[k] = AggInfo[vertex2AggId[k]];
744 nAggregates = NewNAggs;
752 for (
LO i = 0; i < size; i++) {
763 #endif // MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP
#define MueLu_sumAll(rcpComm, in, out)
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.
const RCP< LOVector > & GetProcWinner() const
Returns constant vector that maps local node IDs to owning processor IDs.
#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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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...
const_pointer const_iterator
void AggregateLeftovers(GraphBase const &graph, Aggregates &aggregates) const
Take a partially aggregated graph and complete the aggregation.
virtual const RCP< const Map > GetDomainMap() const =0
const RCP< const Map > GetMap() const
returns (overlapping) map of aggregate/node distribution
virtual Xpetra::global_size_t GetGlobalNumEdges() const =0
Return number of global edges in the graph.
#define MueLu_minAll(rcpComm, in, out)
void SetIsRoot(LO i, bool value=true)
Set root node information.
Teuchos::ArrayRCP< LO > ComputeAggregateSizes(bool forceRecompute=false) const
Compute sizes of aggregates.
void ArbitrateAndCommunicate(Vector &weights, Aggregates &aggregates, const bool perturb) const
This method assigns unknowns to aggregates.
Helper class for providing arbitrated communication across processors.
#define MUELU_UNAGGREGATED
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.
int RemoveSmallAggs(Aggregates &aggregates, int min_size, RCP< Xpetra::Vector< SC, LO, GO, NO > > &distWeights, const MueLu::CoupledAggregationCommHelper< LO, GO, NO > &myWidget) const
Attempt to clean up aggregates that are too small.
virtual const RCP< const Teuchos::Comm< int > > GetComm() const =0
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
MueLu representation of a graph.
LeftoverAggregationAlgorithm()
Constructor.
Timer to be used in non-factories.
bool IsRoot(LO i) const
Returns true if node with given local node id is marked to be a root node.
#define MUELU_PHASE4BUCKETS
#define MUELU_DISTONE_VERTEX_WEIGHT
Exception throws to report errors in the internal logical of the program.
#define MUELU_PENALTYFACTOR
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.