MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_LeftoverAggregationAlgorithm_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP
47#define MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP
48
49#include <Xpetra_Map.hpp>
50#include <Xpetra_MultiVectorFactory.hpp>
51#include <Xpetra_VectorFactory.hpp>
52
54
55#include "MueLu_Aggregates_decl.hpp" // MUELU_UNASSIGNED macro
56#include "MueLu_Utilities_decl.hpp" // MueLu_sumAll macro
57#include "MueLu_GraphBase.hpp"
58#include "MueLu_CoupledAggregationCommHelper.hpp"
59#include "MueLu_Exceptions.hpp"
60#include "MueLu_Monitor.hpp"
61
62#ifdef HAVE_MUELU_TPETRA
63#include <Tpetra_Details_DefaultTypes.hpp>
64#endif
65
66
67namespace MueLu {
68
69 template <class LocalOrdinal, class GlobalOrdinal, class Node>
71 phase3AggCreation_(.5),
72 minNodesPerAggregate_(1)
73 { }
74
75 template <class LocalOrdinal, class GlobalOrdinal, class Node>
77 Monitor m(*this, "AggregateLeftovers");
78
79 my_size_t nVertices = graph.GetNodeNumVertices();
80 int exp_nRows = aggregates.GetMap()->getLocalNumElements(); // Tentative fix... was previously exp_nRows = nVertices + graph.GetNodeNumGhost();
81 int myPid = graph.GetComm()->getRank();
82 my_size_t nAggregates = aggregates.GetNumAggregates();
83
84 int minNodesPerAggregate = GetMinNodesPerAggregate();
85
86 const RCP<const Map> nonUniqueMap = aggregates.GetMap(); //column map of underlying graph
87 const RCP<const Map> uniqueMap = graph.GetDomainMap();
88
89 MueLu::CoupledAggregationCommHelper<LO,GO,NO> myWidget(uniqueMap, nonUniqueMap);
90
91 //TODO JJH We want to skip this call
92 RCP<Xpetra::Vector<SC,LO,GO,NO> > distWeights = Xpetra::VectorFactory<SC,LO,GO,NO>::Build(nonUniqueMap);
93
94 // Aggregated vertices not "definitively" assigned to processors are
95 // arbitrated by ArbitrateAndCommunicate(). There is some
96 // additional logic to prevent losing root nodes in arbitration.
97 {
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);
101
102 for (size_t i=0;i<nonUniqueMap->getLocalNumElements();i++) {
103 if (procWinner[i] == MUELU_UNASSIGNED) {
104 if (vertex2AggId[i] != MUELU_UNAGGREGATED) {
105 weights[i] = 1.;
106 if (aggregates.IsRoot(i)) weights[i] = 2.;
107 }
108 }
109 }
110
111 // views on distributed vectors are freed here.
112 }
113
114 //TODO JJH We want to skip this call
115 myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
116 // All tentatively assigned vertices are now definitive
117
118 // Tentatively assign any vertex (ghost or local) which neighbors a root
119 // to the aggregate associated with the root.
120 {
121 ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
122 ArrayRCP<const LO> procWinner = aggregates.GetProcWinner()->getData(0);
123 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
124
125 for (my_size_t i = 0; i < nVertices; i++) {
126 if ( aggregates.IsRoot(i) && (procWinner[i] == myPid) ) {
127
128 // neighOfINode is the neighbor node list of node 'i'.
129 ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
130
131 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
132 int colj = *it;
133 if (vertex2AggId[colj] == MUELU_UNAGGREGATED) {
134 weights[colj]= 1.;
135 vertex2AggId[colj] = vertex2AggId[i];
136 }
137 }
138 }
139 }
140
141 // views on distributed vectors are freed here.
142 }
143
144 //TODO JJH We want to skip this call
145 myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
146 // All tentatively assigned vertices are now definitive
147
148 // Record the number of aggregated vertices
149 GO total_phase_one_aggregated = 0;
150 {
151 ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
152
153 GO phase_one_aggregated = 0;
154 for (my_size_t i = 0; i < nVertices; i++) {
155 if (vertex2AggId[i] != MUELU_UNAGGREGATED)
156 phase_one_aggregated++;
157 }
158
159 MueLu_sumAll(graph.GetComm(), phase_one_aggregated, total_phase_one_aggregated);
160
161 GO local_nVertices = nVertices, total_nVertices = 0;
162 MueLu_sumAll(graph.GetComm(), local_nVertices, total_nVertices);
163
164 /* Among unaggregated points, see if we can make a reasonable size */
165 /* aggregate out of it. We do this by looking at neighbors and seeing */
166 /* how many are unaggregated and on my processor. Loosely, */
167 /* base the number of new aggregates created on the percentage of */
168 /* unaggregated nodes. */
169
170 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
171
172 double factor = 1.;
173 factor = ((double) total_phase_one_aggregated)/((double)(total_nVertices + 1));
174 factor = pow(factor, GetPhase3AggCreation());
175
176 for (my_size_t i = 0; i < nVertices; i++) {
177 if (vertex2AggId[i] == MUELU_UNAGGREGATED)
178 {
179
180 // neighOfINode is the neighbor node list of node 'iNode'.
181 ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
182 int rowi_N = neighOfINode.size();
183
184 int nonaggd_neighbors = 0;
185 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
186 int colj = *it;
187 if (vertex2AggId[colj] == MUELU_UNAGGREGATED && colj < nVertices)
188 nonaggd_neighbors++;
189 }
190 if ( (nonaggd_neighbors > minNodesPerAggregate) &&
191 (((double) nonaggd_neighbors)/((double) rowi_N) > factor))
192 {
193 vertex2AggId[i] = (nAggregates)++;
194 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
195 int colj = *it;
196 if (vertex2AggId[colj]==MUELU_UNAGGREGATED) {
197 vertex2AggId[colj] = vertex2AggId[i];
198 if (colj < nVertices) weights[colj] = 2.;
199 else weights[colj] = 1.;
200 }
201 }
202 aggregates.SetIsRoot(i);
203 weights[i] = 2.;
204 }
205 }
206 } // for (i = 0; i < nVertices; i++)
207
208 // views on distributed vectors are freed here.
209 }
210
211 //TODO JJH We want to skip this call
212 myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
213 //All tentatively assigned vertices are now definitive
214
215 if (IsPrint(Statistics1)) {
216 GO Nphase1_agg = nAggregates;
217 GO total_aggs;
218
219 MueLu_sumAll(graph.GetComm(), Nphase1_agg, total_aggs);
220
221 GetOStream(Statistics1) << "Phase 1 - nodes aggregated = " << total_phase_one_aggregated << std::endl;
222 GetOStream(Statistics1) << "Phase 1 - total aggregates = " << total_aggs << std::endl;
223
224 GO i = nAggregates - Nphase1_agg;
225 { GO ii; MueLu_sumAll(graph.GetComm(),i,ii); i = ii; }
226 GetOStream(Statistics1) << "Phase 3 - additional aggregates = " << i << std::endl;
227 }
228
229 // Determine vertices that are not shared by setting Temp to all ones
230 // and doing NonUnique2NonUnique(..., ADD). This sums values of all
231 // local copies associated with each Gid. Thus, sums > 1 are shared.
232
233 // std::cout << "exp_nrows=" << exp_nRows << " (nVertices= " << nVertices << ", numGhost=" << graph.GetNodeNumGhost() << ")" << std::endl;
234 // std::cout << "nonUniqueMap=" << nonUniqueMap->getLocalNumElements() << std::endl;
235
236 RCP<Xpetra::Vector<SC,LO,GO,NO> > temp_ = Xpetra::VectorFactory<SC,LO,GO,NO> ::Build(nonUniqueMap,false); //no need to zero out vector in ctor
237 temp_->putScalar(1.);
238
239 RCP<Xpetra::Vector<SC,LO,GO,NO> > tempOutput_ = Xpetra::VectorFactory<SC,LO,GO,NO> ::Build(nonUniqueMap);
240
241 myWidget.NonUnique2NonUnique(*temp_, *tempOutput_, Xpetra::ADD);
242
243 std::vector<bool> gidNotShared(exp_nRows);
244 {
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;
249 else
250 gidNotShared[i] = true;
251 }
252 }
253
254 // Phase 4.
255 double nAggregatesTarget;
256 nAggregatesTarget = ((double) uniqueMap->getGlobalNumElements())* (((double) uniqueMap->getGlobalNumElements())/ ((double) graph.GetGlobalNumEdges()));
257
258 GO nAggregatesLocal=nAggregates, nAggregatesGlobal; MueLu_sumAll(graph.GetComm(), nAggregatesLocal, nAggregatesGlobal);
259
260 LO minNAggs; MueLu_minAll(graph.GetComm(), nAggregates, minNAggs);
261 LO maxNAggs; MueLu_maxAll(graph.GetComm(), nAggregates, maxNAggs);
262
263 //
264 // Only do this phase if things look really bad. THIS
265 // CODE IS PRETTY EXPERIMENTAL
266 //
267#define MUELU_PHASE4BUCKETS 6
268 if ((nAggregatesGlobal < graph.GetComm()->getSize()) &&
269 (2.5*nAggregatesGlobal < nAggregatesTarget) &&
270 (minNAggs ==0) && (maxNAggs <= 1)) {
271
272 // Modify seed of the random algorithm used by temp_->randomize()
273 {
274 typedef Teuchos::ScalarTraits<double> scalarTrait; // temp_ is of type double.
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()));
279 }
280
281 temp_->randomize();
282
283 ArrayRCP<SC> temp = temp_->getDataNonConst(0);
284
285 // build a list of candidate root nodes (vertices not adjacent
286 // to aggregated vertices)
287
288 my_size_t nCandidates = 0;
289 global_size_t nCandidatesGlobal;
290
291 ArrayRCP<LO> candidates = Teuchos::arcp<LO>(nVertices+1);
292
293 double priorThreshold = 0.;
294 for (int kkk = 0; kkk < MUELU_PHASE4BUCKETS; kkk++) {
295
296 {
297 ArrayRCP<const LO> vertex2AggId = aggregates.GetVertex2AggId()->getData(0);
298 ArrayView<const LO> vertex2AggIdView = vertex2AggId();
299 RootCandidates(nVertices, vertex2AggIdView, graph, candidates, nCandidates, nCandidatesGlobal);
300 // views on distributed vectors are freed here.
301 }
302
303 double nTargetNewGuys = nAggregatesTarget - nAggregatesGlobal;
304 double threshold = priorThreshold + (1. - priorThreshold)*nTargetNewGuys/(nCandidatesGlobal + .001);
305
306 threshold = (threshold*(kkk+1.))/((double) MUELU_PHASE4BUCKETS);
307 priorThreshold = threshold;
308
309 {
310 ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
311 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
312
313 for (int k = 0; k < nCandidates; k++ ) {
314 int i = candidates[k];
315 if ((vertex2AggId[i] == MUELU_UNAGGREGATED) && (fabs(temp[i]) < threshold)) {
316 // Note: priorThreshold <= fabs(temp[i]) <= 1
317
318 // neighOfINode is the neighbor node list of node 'iNode'.
319 ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
320
321 if (neighOfINode.size() > minNodesPerAggregate) { //TODO: check if this test is exactly was we want to do
322 int count = 0;
323 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
324 int Adjacent = *it;
325 // This might not be true if someone close to i
326 // is chosen as a root via fabs(temp[]) < Threshold
327 if (vertex2AggId[Adjacent] == MUELU_UNAGGREGATED){
328 count++;
329 vertex2AggId[Adjacent] = nAggregates;
330 weights[Adjacent] = 1.;
331 }
332 }
333 if (count >= minNodesPerAggregate) {
334 vertex2AggId[i] = nAggregates++;
335 weights[i] = 2.;
336 aggregates.SetIsRoot(i);
337 }
338 else { // undo things
339 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
340 int Adjacent = *it;
341 if (vertex2AggId[Adjacent] == nAggregates){
342 vertex2AggId[Adjacent] = MUELU_UNAGGREGATED;
343 weights[Adjacent] = 0.;
344 }
345 }
346 }
347 }
348 }
349 }
350 // views on distributed vectors are freed here.
351 }
352 //TODO JJH We want to skip this call
353 myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
354 // All tentatively assigned vertices are now definitive
355 nAggregatesLocal=nAggregates;
356 MueLu_sumAll(graph.GetComm(), nAggregatesLocal, nAggregatesGlobal);
357
358 // check that there are no aggregates sizes below minNodesPerAggregate
359
360 aggregates.SetNumAggregates(nAggregates);
361
362 RemoveSmallAggs(aggregates, minNodesPerAggregate, distWeights, myWidget);
363
364 nAggregates = aggregates.GetNumAggregates();
365 } // one possibility
366 }
367
368 // Initialize things for Phase 5. This includes building the transpose
369 // of the matrix ONLY for transposed rows that correspond to unaggregted
370 // ghost vertices. Further, the transpose is only a local transpose.
371 // Nonzero edges which exist on other processors are not represented.
372
373
374 int observedNAgg=-1; //number of aggregates that contain vertices on this process
375
376 {
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];
382 observedNAgg++;
383 }
384
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);
388
389 for (int i = 0; i < exp_nRows; i++) Mark[i] = MUELU_DISTONE_VERTEX_WEIGHT;
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;
392
393 // Grab the transpose matrix graph for unaggregated ghost vertices.
394 // a) count the number of nonzeros per row in the transpose
395 std::vector<int> RowPtr(exp_nRows+1-nVertices);
396 //{
397 ArrayRCP<const LO> vertex2AggIdCst = aggregates.GetVertex2AggId()->getData(0);
398
399 for (int i = nVertices; i < exp_nRows; i++) RowPtr[i-nVertices] = 0;
400 for (int i = 0; i < nVertices; i++) {
401
402 // neighOfINode is the neighbor node list of node 'iNode'.
403 ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
404
405 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
406 int j = *it;
407 if ( (j >= nVertices) && (vertex2AggIdCst[j] == MUELU_UNAGGREGATED)){
408 RowPtr[j-nVertices]++;
409 }
410 }
411 }
412
413 // b) Convert RowPtr[i] to point to 1st first nnz spot in row i.
414
415 int iSum = 0, iTemp;
416 for (int i = nVertices; i < exp_nRows; i++) {
417 iTemp = RowPtr[i-nVertices];
418 RowPtr[i-nVertices] = iSum;
419 iSum += iTemp;
420 }
421 RowPtr[exp_nRows-nVertices] = iSum;
422 std::vector<LO> cols(iSum+1);
423
424 // c) Traverse matrix and insert entries in proper location.
425 for (int i = 0; i < nVertices; i++) {
426
427 // neighOfINode is the neighbor node list of node 'iNode'.
428 ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
429
430 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
431 int j = *it;
432 if ( (j >= nVertices) && (vertex2AggIdCst[j] == MUELU_UNAGGREGATED)){
433 cols[RowPtr[j-nVertices]++] = i;
434 }
435 }
436 }
437
438 // d) RowPtr[i] points to beginning of row i+1 so shift by one location.
439 for (int i = exp_nRows; i > nVertices; i--)
440 RowPtr[i-nVertices] = RowPtr[i-1-nVertices];
441 RowPtr[0] = 0;
442
443 // views on distributed vectors are freed here.
444 vertex2AggIdCst = Teuchos::null;
445 //}
446
447 int bestScoreCutoff;
448 int thresholds[10] = {300,200,100,50,25,13,7,4,2,0};
449
450 // Stick unaggregated vertices into existing aggregates as described above.
451
452 {
453 int ncalls=0;
454
455 for (int kk = 0; kk < 10; kk += 2) {
456 bestScoreCutoff = thresholds[kk];
457
458 ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
459 ArrayRCP<const LO> procWinner = aggregates.GetProcWinner()->getData(0);
460 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
461
462 for (int i = 0; i < exp_nRows; i++) {
463
464 if (vertex2AggId[i] == MUELU_UNAGGREGATED) {
465
466 // neighOfINode is the neighbor node list of node 'iNode'.
467 ArrayView<const LO> neighOfINode;
468
469 // Grab neighboring vertices which is either in graph for local ids
470 // or sits in transposed fragment just constructed above for ghosts.
471 if (i < nVertices) {
472 neighOfINode = graph.getNeighborVertices(i);
473 }
474 else {
475 LO *rowi_col = NULL, rowi_N;
476 rowi_col = &(cols[RowPtr[i-nVertices]]);
477 rowi_N = RowPtr[i+1-nVertices] - RowPtr[i-nVertices];
478
479 neighOfINode = ArrayView<const LO>(rowi_col, rowi_N);
480 }
481 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
482 int Adjacent = *it;
483 int AdjacentAgg = vertex2AggId[Adjacent];
484
485 //Adjacent is aggregated and either I own the aggregate
486 // or I could own the aggregate after arbitration.
487 if ((AdjacentAgg != MUELU_UNAGGREGATED) &&
488 ((procWinner[Adjacent] == myPid) ||
489 (procWinner[Adjacent] == MUELU_UNASSIGNED))){
490 SumOfMarks[AdjacentAgg] += Mark[Adjacent];
491 }
492 }
493 int best_score = MUELU_NOSCORE;
494 int best_agg = -1;
495 int BestMark = -1;
496 bool cannotLoseAllFriends=false; // Used to address possible loss of vertices in arbitration of shared nodes discussed above. (Initialized to false only to avoid a compiler warning).
497
498 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
499 int Adjacent = *it;
500 int AdjacentAgg = vertex2AggId[Adjacent];
501 //Adjacent is unaggregated, has some value and no
502 //other processor has definitively claimed him
503 if ((AdjacentAgg != MUELU_UNAGGREGATED) &&
504 (SumOfMarks[AdjacentAgg] != 0) &&
505 ((procWinner[Adjacent] == myPid) ||
506 (procWinner[Adjacent] == MUELU_UNASSIGNED ))) {
507
508 // first figure out the penalty associated with
509 // AdjacentAgg having already been incremented
510 // during this phase, then compute score.
511
512 double penalty = (double) (INCR_SCALING*agg_incremented[AdjacentAgg]);
513 if (penalty > MUELU_PENALTYFACTOR*((double)SumOfMarks[AdjacentAgg]))
514 penalty = MUELU_PENALTYFACTOR*((double)SumOfMarks[AdjacentAgg]);
515 int score = SumOfMarks[AdjacentAgg]- ((int) floor(penalty));
516
517 if (score > best_score) {
518 best_agg = AdjacentAgg;
519 best_score = score;
520 BestMark = Mark[Adjacent];
521 cannotLoseAllFriends = false;
522
523 // This address issue mentioned above by checking whether
524 // Adjacent could be lost in arbitration. weight==0 means that
525 // Adjacent was not set during this loop of Phase 5 (and so it
526 // has already undergone arbitration). GidNotShared == true
527 // obviously implies that Adjacent cannot be lost to arbitration
528 if ((weights[Adjacent]== 0.) || (gidNotShared[Adjacent] == true))
529 cannotLoseAllFriends = true;
530 }
531 // Another vertex within current best aggregate found.
532 // We should have (best_score == score). We need to see
533 // if we can improve BestMark and cannotLoseAllFriends.
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];
538 }
539 }
540 }
541 // Clean up
542 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
543 int Adjacent = *it;
544 int AdjacentAgg = vertex2AggId[Adjacent];
545 if (AdjacentAgg >= 0) SumOfMarks[AdjacentAgg] = 0;
546 }
547 // Tentatively assign vertex to best_agg.
548 if ( (best_score >= bestScoreCutoff) && (cannotLoseAllFriends)) {
549
550 TEUCHOS_TEST_FOR_EXCEPTION(best_agg == -1 || BestMark == -1, MueLu::Exceptions::RuntimeError, "MueLu::CoupledAggregationFactory internal error"); // should never happen
551
552 vertex2AggId[i] = best_agg;
553 weights[i] = best_score;
554 agg_incremented[best_agg]++;
555 Mark[i] = (int) ceil( ((double) BestMark)/2.);
556 }
557 }
558
559 // views on distributed vectors are freed here.
560 }
561
562 vertex2AggId = Teuchos::null;
563 procWinner = Teuchos::null;
564 weights = Teuchos::null;
565
566 ++ncalls;
567 //TODO JJH We want to skip this call
568 myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
569 // All tentatively assigned vertices are now definitive
570 }
571
572 // if (graph.GetComm()->getRank()==0)
573 // std::cout << "#calls to Arb&Comm=" << ncalls << std::endl;
574 }
575
576 // Phase 6: Aggregate remain unaggregated vertices and try at all costs
577 // to avoid small aggregates.
578 // One case where we can find ourselves in this situation
579 // is if all vertices vk adjacent to v have already been
580 // put in other processor's aggregates and v does not have
581 // a direct connection to a local vertex in any of these
582 // aggregates.
583
584 int Nleftover = 0, Nsingle = 0;
585 {
586
587 ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
588 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
589 ArrayRCP<const LO> procWinner = aggregates.GetProcWinner()->getData(0);
590
591 int count = 0;
592 for (my_size_t i = 0; i < nVertices; i++) {
593 if (vertex2AggId[i] == MUELU_UNAGGREGATED) {
594 Nleftover++;
595
596 // neighOfINode is the neighbor node list of node 'iNode'.
597 ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
598
599 // We don't want too small of an aggregate. So lets see if there is an
600 // unaggregated neighbor that we can also put with this vertex
601
602 vertex2AggId[i] = nAggregates;
603 weights[i] = 1.;
604 if (count == 0) aggregates.SetIsRoot(i);
605 count++;
606 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
607 int j = *it;
608 if ((j != i)&&(vertex2AggId[j] == MUELU_UNAGGREGATED)&&
609 (j < nVertices)) {
610 vertex2AggId[j] = nAggregates;
611 weights[j] = 1.;
612 count++;
613 }
614 }
615 if ( count >= minNodesPerAggregate) {
616 nAggregates++;
617 count = 0;
618 }
619 }
620 }
621
622 // We have something which is under minNodesPerAggregate when
623 if (count != 0) {
624#ifdef FIXME
625 // Can stick small aggregate with 0th aggregate?
626 if (nAggregates > 0) {
627 for (my_size_t i = 0; i < nVertices; i++) {
628 if ((vertex2AggId[i] == nAggregates) && (procWinner[i] == myPid)) {
629 vertex2AggId[i] = 0;
630 aggregates.SetIsRoot(i,false);
631 }
632 }
633 }
634 else {
635 Nsingle++;
636 nAggregates++;
637 }
638#else
639 // Can stick small aggregate with 0th aggregate?
640 if (nAggregates > 0) {
641 for (my_size_t i = 0; i < nVertices; i++) {
642 // TW: This is not a real fix. This may produce ugly bad aggregates!
643 // I removed the procWinner[i] == myPid check. it makes no sense to me since
644 // it leaves vertex2AggId[i] == nAggregates -> crash in ComputeAggregateSizes().
645 // Maybe it's better to add the leftovers to the last generated agg on the current proc.
646 // The best solution would be to add them to the "next"/nearest aggregate, that may be
647 // on an other processor
648 if (vertex2AggId[i] == nAggregates) {
649 vertex2AggId[i] = nAggregates-1; //0;
650 aggregates.SetIsRoot(i,false);
651 }
652 }
653 }
654 else {
655 Nsingle++;
656 nAggregates++;
657 }
658#endif
659 }
660
661 // views on distributed vectors are freed here.
662 }
663
664 //TODO JJH We want to skip this call
665 myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, false);
666
667 if (IsPrint(Statistics1)) {
668 GO total_Nsingle=0; MueLu_sumAll(graph.GetComm(), (GO)Nsingle, total_Nsingle);
669 GO total_Nleftover=0; MueLu_sumAll(graph.GetComm(), (GO)Nleftover, total_Nleftover);
670 // GO total_aggs; MueLu_sumAll(graph.GetComm(), (GO)nAggregates, total_aggs);
671 // GetOStream(Statistics1) << "Phase 6 - total aggregates = " << total_aggs << std::endl;
672 GetOStream(Statistics1) << "Phase 6 - leftovers = " << total_Nleftover << " and singletons = " << total_Nsingle << std::endl;
673 }
674
675 aggregates.SetNumAggregates(nAggregates);
676
677 } //AggregateLeftovers
678
679 template <class LocalOrdinal, class GlobalOrdinal, class Node>
681 ArrayView<const LO> & vertex2AggId, GraphBase const &graph,
682 ArrayRCP<LO> &candidates, my_size_t &nCandidates, global_size_t &nCandidatesGlobal) const
683 {
684 nCandidates = 0;
685
686 for (my_size_t i = 0; i < nVertices; i++ ) {
687 if (vertex2AggId[i] == MUELU_UNAGGREGATED) {
688 bool noAggdNeighbors = true;
689
690 // neighOfINode is the neighbor node list of node 'iNode'.
691 ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
692
693 for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
694 int adjacent = *it;
695 if (vertex2AggId[adjacent] != MUELU_UNAGGREGATED)
696 noAggdNeighbors = false;
697 }
698 if (noAggdNeighbors == true) candidates[nCandidates++] = i;
699 }
700 }
701
702 MueLu_sumAll(graph.GetComm(), (GO)nCandidates, nCandidatesGlobal);
703
704 } //RootCandidates
705
706 template <class LocalOrdinal, class GlobalOrdinal, class Node>
708 RCP<Xpetra::Vector<SC,LO,GO,NO> > & distWeights, const MueLu::CoupledAggregationCommHelper<LO,GO,NO> & myWidget) const {
709 int myPid = aggregates.GetMap()->getComm()->getRank();
710
711 LO nAggregates = aggregates.GetNumAggregates();
712
713 ArrayRCP<LO> procWinner = aggregates.GetProcWinner()->getDataNonConst(0);
714 ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
715 LO size = procWinner.size();
716
717 //ArrayRCP<int> AggInfo = Teuchos::arcp<int>(nAggregates+1);
718 //aggregates.ComputeAggSizes(AggInfo);
719 ArrayRCP<LO> AggInfo = aggregates.ComputeAggregateSizes();
720
721 ArrayRCP<SC> weights = distWeights->getDataNonConst(0);
722
723 // Make a list of all aggregates indicating New AggId
724 // Use AggInfo array for this.
725
726 LO NewNAggs = 0;
727 for (LO i = 0; i < nAggregates; i++) {
728 if ( AggInfo[i] < min_size) {
729 AggInfo[i] = MUELU_UNAGGREGATED;
730 }
731 else AggInfo[i] = NewNAggs++;
732 }
733
734 for (LO k = 0; k < size; k++ ) {
735 if (procWinner[k] == myPid) {
736 if (vertex2AggId[k] != MUELU_UNAGGREGATED) {
737 vertex2AggId[k] = AggInfo[vertex2AggId[k]];
738 weights[k] = 1.;
739 }
740 if (vertex2AggId[k] == MUELU_UNAGGREGATED)
741 aggregates.SetIsRoot(k,false);
742 }
743 }
744 nAggregates = NewNAggs;
745
746 //TODO JJH We want to skip this call
747 myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
748 // All tentatively assigned vertices are now definitive
749
750 // procWinner is not set correctly for aggregates which have
751 // been eliminated
752 for (LO i = 0; i < size; i++) {
753 if (vertex2AggId[i] == MUELU_UNAGGREGATED)
754 procWinner[i] = MUELU_UNASSIGNED;
755 }
756 aggregates.SetNumAggregates(nAggregates);
757
758 return 0; //TODO
759 } //RemoveSmallAggs
760
761} //namespace MueLu
762
763#endif // MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP
#define MUELU_UNAGGREGATED
#define MUELU_UNASSIGNED
#define MUELU_PHASE4BUCKETS
#define MueLu_maxAll(rcpComm, in, out)
#define MueLu_sumAll(rcpComm, in, out)
#define MueLu_minAll(rcpComm, in, out)
Container class for aggregation information.
Teuchos::ArrayRCP< LO > ComputeAggregateSizes(bool forceRecompute=false) const
Compute sizes of aggregates.
const RCP< const Map > GetMap() const
returns (overlapping) map of aggregate/node distribution
void SetIsRoot(LO i, bool value=true)
Set root node information.
bool IsRoot(LO i) const
Returns true if node with given local node id is marked to be a root node.
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
const RCP< LOVector > & GetProcWinner() const
Returns constant vector that maps local node IDs to owning processor IDs.
LO GetNumAggregates() const
returns the number of aggregates of the current processor. Note: could/should be renamed to GetNumLoc...
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.
Helper class for providing arbitrated communication across processors.
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...
void ArbitrateAndCommunicate(Vector &weights, Aggregates &aggregates, const bool perturb) const
This method assigns unknowns to aggregates.
Exception throws to report errors in the internal logical of the program.
MueLu representation of a graph.
virtual const RCP< const Teuchos::Comm< int > > GetComm() const =0
virtual const RCP< const Map > GetDomainMap() const =0
virtual Xpetra::global_size_t GetGlobalNumEdges() const =0
Return number of global edges in the graph.
virtual Teuchos::ArrayView< const LocalOrdinal > getNeighborVertices(LocalOrdinal v) const =0
Return the list of vertices adjacent to the vertex 'v'.
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
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.
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.
void AggregateLeftovers(GraphBase const &graph, Aggregates &aggregates) const
Take a partially aggregated graph and complete the aggregation.
Timer to be used in non-factories.
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.