77 Monitor m(*
this,
"AggregateLeftovers");
80 int exp_nRows = aggregates.
GetMap()->getLocalNumElements();
81 int myPid = graph.
GetComm()->getRank();
84 int minNodesPerAggregate = GetMinNodesPerAggregate();
86 const RCP<const Map> nonUniqueMap = aggregates.
GetMap();
92 RCP<Xpetra::Vector<SC,LO,GO,NO> > distWeights = Xpetra::VectorFactory<SC,LO,GO,NO>::Build(nonUniqueMap);
98 ArrayRCP<const LO> vertex2AggId = aggregates.
GetVertex2AggId()->getData(0);
99 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
100 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
102 for (
size_t i=0;i<nonUniqueMap->getLocalNumElements();i++) {
106 if (aggregates.
IsRoot(i)) weights[i] = 2.;
121 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
122 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
123 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
125 for (
my_size_t i = 0; i < nVertices; i++) {
126 if ( aggregates.
IsRoot(i) && (procWinner[i] == myPid) ) {
131 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
135 vertex2AggId[colj] = vertex2AggId[i];
149 GO total_phase_one_aggregated = 0;
151 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(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;
170 ArrayRCP<SC> weights = distWeights->getDataNonConst(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;
185 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
190 if ( (nonaggd_neighbors > minNodesPerAggregate) &&
191 (((double) nonaggd_neighbors)/((double) rowi_N) > factor))
193 vertex2AggId[i] = (nAggregates)++;
194 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
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;
236 RCP<Xpetra::Vector<SC,LO,GO,NO> > temp_ = Xpetra::VectorFactory<SC,LO,GO,NO> ::Build(nonUniqueMap,
false);
237 temp_->putScalar(1.);
239 RCP<Xpetra::Vector<SC,LO,GO,NO> > tempOutput_ = Xpetra::VectorFactory<SC,LO,GO,NO> ::Build(nonUniqueMap);
243 std::vector<bool> gidNotShared(exp_nRows);
245 ArrayRCP<const SC> tempOutput = tempOutput_->getData(0);
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)) {
274 typedef Teuchos::ScalarTraits<double> scalarTrait;
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()));
283 ArrayRCP<SC> temp = temp_->getDataNonConst(0);
291 ArrayRCP<LO> candidates = Teuchos::arcp<LO>(nVertices+1);
293 double priorThreshold = 0.;
297 ArrayRCP<const LO> vertex2AggId = aggregates.
GetVertex2AggId()->getData(0);
298 ArrayView<const LO> vertex2AggIdView = vertex2AggId();
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;
310 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
311 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
313 for (
int k = 0; k < nCandidates; k++ ) {
314 int i = candidates[k];
321 if (neighOfINode.size() > minNodesPerAggregate) {
323 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
329 vertex2AggId[Adjacent] = nAggregates;
330 weights[Adjacent] = 1.;
333 if (count >= minNodesPerAggregate) {
334 vertex2AggId[i] = nAggregates++;
339 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
341 if (vertex2AggId[Adjacent] == nAggregates){
343 weights[Adjacent] = 0.;
355 nAggregatesLocal=nAggregates;
362 RemoveSmallAggs(aggregates, minNodesPerAggregate, distWeights, myWidget);
377 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
378 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
379 for(LO k = 0; k < vertex2AggId.size(); ++k )
380 if(vertex2AggId[k]>observedNAgg)
381 observedNAgg=vertex2AggId[k];
385 ArrayRCP<int> Mark = Teuchos::arcp<int>(exp_nRows+1);
386 ArrayRCP<int> agg_incremented = Teuchos::arcp<int>(observedNAgg);
387 ArrayRCP<int> SumOfMarks = 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);
397 ArrayRCP<const LO> vertex2AggIdCst = aggregates.
GetVertex2AggId()->getData(0);
399 for (
int i = nVertices; i < exp_nRows; i++) RowPtr[i-nVertices] = 0;
400 for (
int i = 0; i < nVertices; i++) {
405 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
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++) {
430 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
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];
458 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
459 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
460 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
462 for (
int i = 0; i < exp_nRows; i++) {
467 ArrayView<const LO> neighOfINode;
475 LO *rowi_col = NULL, rowi_N;
476 rowi_col = &(cols[RowPtr[i-nVertices]]);
477 rowi_N = RowPtr[i+1-nVertices] - RowPtr[i-nVertices];
479 neighOfINode = ArrayView<const LO>(rowi_col, rowi_N);
481 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
483 int AdjacentAgg = vertex2AggId[Adjacent];
488 ((procWinner[Adjacent] == myPid) ||
490 SumOfMarks[AdjacentAgg] += Mark[Adjacent];
496 bool cannotLoseAllFriends=
false;
498 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
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];
542 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
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;
587 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
588 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
589 ArrayRCP<const LO> procWinner = aggregates.
GetProcWinner()->getData(0);
592 for (
my_size_t i = 0; i < nVertices; i++) {
602 vertex2AggId[i] = nAggregates;
606 for (
typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
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;