MueLu  Version of the Day
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_VectorFactory.hpp>
51 
53 
54 #include "MueLu_Aggregates_decl.hpp" // MUELU_UNASSIGNED macro
55 #include "MueLu_Utilities_decl.hpp" // MueLu_sumAll macro
56 #include "MueLu_GraphBase.hpp"
57 #include "MueLu_CoupledAggregationCommHelper.hpp"
58 #include "MueLu_Exceptions.hpp"
59 #include "MueLu_Monitor.hpp"
60 
61 namespace MueLu {
62 
63  template <class LocalOrdinal, class GlobalOrdinal, class Node>
65  phase3AggCreation_(.5),
66  minNodesPerAggregate_(1)
67  { }
68 
69  template <class LocalOrdinal, class GlobalOrdinal, class Node>
71  Monitor m(*this, "AggregateLeftovers");
72 
73  my_size_t nVertices = graph.GetNodeNumVertices();
74  int exp_nRows = aggregates.GetMap()->getNodeNumElements(); // Tentative fix... was previously exp_nRows = nVertices + graph.GetNodeNumGhost();
75  int myPid = graph.GetComm()->getRank();
76  my_size_t nAggregates = aggregates.GetNumAggregates();
77 
78  int minNodesPerAggregate = GetMinNodesPerAggregate();
79 
80  const RCP<const Map> nonUniqueMap = aggregates.GetMap(); //column map of underlying graph
81  const RCP<const Map> uniqueMap = graph.GetDomainMap();
82 
83  MueLu::CoupledAggregationCommHelper<LO,GO,NO> myWidget(uniqueMap, nonUniqueMap);
84 
85  //TODO JJH We want to skip this call
86  RCP<Xpetra::Vector<double,LO,GO,NO> > distWeights = Xpetra::VectorFactory<double,LO,GO,NO>::Build(nonUniqueMap);
87 
88  // Aggregated vertices not "definitively" assigned to processors are
89  // arbitrated by ArbitrateAndCommunicate(). There is some
90  // additional logic to prevent losing root nodes in arbitration.
91  {
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);
95 
96  for (size_t i=0;i<nonUniqueMap->getNodeNumElements();i++) {
97  if (procWinner[i] == MUELU_UNASSIGNED) {
98  if (vertex2AggId[i] != MUELU_UNAGGREGATED) {
99  weights[i] = 1.;
100  if (aggregates.IsRoot(i)) weights[i] = 2.;
101  }
102  }
103  }
104 
105  // views on distributed vectors are freed here.
106  }
107 
108  //TODO JJH We want to skip this call
109  myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
110  // All tentatively assigned vertices are now definitive
111 
112  // Tentatively assign any vertex (ghost or local) which neighbors a root
113  // to the aggregate associated with the root.
114  {
115  ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
116  ArrayRCP<const LO> procWinner = aggregates.GetProcWinner()->getData(0);
117  ArrayRCP<double> weights = distWeights->getDataNonConst(0);
118 
119  for (my_size_t i = 0; i < nVertices; i++) {
120  if ( aggregates.IsRoot(i) && (procWinner[i] == myPid) ) {
121 
122  // neighOfINode is the neighbor node list of node 'i'.
123  ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
124 
125  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
126  int colj = *it;
127  if (vertex2AggId[colj] == MUELU_UNAGGREGATED) {
128  weights[colj]= 1.;
129  vertex2AggId[colj] = vertex2AggId[i];
130  }
131  }
132  }
133  }
134 
135  // views on distributed vectors are freed here.
136  }
137 
138  //TODO JJH We want to skip this call
139  myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
140  // All tentatively assigned vertices are now definitive
141 
142  // Record the number of aggregated vertices
143  GO total_phase_one_aggregated = 0;
144  {
145  ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
146 
147  GO phase_one_aggregated = 0;
148  for (my_size_t i = 0; i < nVertices; i++) {
149  if (vertex2AggId[i] != MUELU_UNAGGREGATED)
150  phase_one_aggregated++;
151  }
152 
153  MueLu_sumAll(graph.GetComm(), phase_one_aggregated, total_phase_one_aggregated);
154 
155  GO local_nVertices = nVertices, total_nVertices = 0;
156  MueLu_sumAll(graph.GetComm(), local_nVertices, total_nVertices);
157 
158  /* Among unaggregated points, see if we can make a reasonable size */
159  /* aggregate out of it. We do this by looking at neighbors and seeing */
160  /* how many are unaggregated and on my processor. Loosely, */
161  /* base the number of new aggregates created on the percentage of */
162  /* unaggregated nodes. */
163 
164  ArrayRCP<double> weights = distWeights->getDataNonConst(0);
165 
166  double factor = 1.;
167  factor = ((double) total_phase_one_aggregated)/((double)(total_nVertices + 1));
168  factor = pow(factor, GetPhase3AggCreation());
169 
170  for (my_size_t i = 0; i < nVertices; i++) {
171  if (vertex2AggId[i] == MUELU_UNAGGREGATED)
172  {
173 
174  // neighOfINode is the neighbor node list of node 'iNode'.
175  ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
176  int rowi_N = neighOfINode.size();
177 
178  int nonaggd_neighbors = 0;
179  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
180  int colj = *it;
181  if (vertex2AggId[colj] == MUELU_UNAGGREGATED && colj < nVertices)
182  nonaggd_neighbors++;
183  }
184  if ( (nonaggd_neighbors > minNodesPerAggregate) &&
185  (((double) nonaggd_neighbors)/((double) rowi_N) > factor))
186  {
187  vertex2AggId[i] = (nAggregates)++;
188  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
189  int colj = *it;
190  if (vertex2AggId[colj]==MUELU_UNAGGREGATED) {
191  vertex2AggId[colj] = vertex2AggId[i];
192  if (colj < nVertices) weights[colj] = 2.;
193  else weights[colj] = 1.;
194  }
195  }
196  aggregates.SetIsRoot(i);
197  weights[i] = 2.;
198  }
199  }
200  } // for (i = 0; i < nVertices; i++)
201 
202  // views on distributed vectors are freed here.
203  }
204 
205  //TODO JJH We want to skip this call
206  myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
207  //All tentatively assigned vertices are now definitive
208 
209  if (IsPrint(Statistics1)) {
210  GO Nphase1_agg = nAggregates;
211  GO total_aggs;
212 
213  MueLu_sumAll(graph.GetComm(), Nphase1_agg, total_aggs);
214 
215  GetOStream(Statistics1) << "Phase 1 - nodes aggregated = " << total_phase_one_aggregated << std::endl;
216  GetOStream(Statistics1) << "Phase 1 - total aggregates = " << total_aggs << std::endl;
217 
218  GO i = nAggregates - Nphase1_agg;
219  { GO ii; MueLu_sumAll(graph.GetComm(),i,ii); i = ii; }
220  GetOStream(Statistics1) << "Phase 3 - additional aggregates = " << i << std::endl;
221  }
222 
223  // Determine vertices that are not shared by setting Temp to all ones
224  // and doing NonUnique2NonUnique(..., ADD). This sums values of all
225  // local copies associated with each Gid. Thus, sums > 1 are shared.
226 
227  // std::cout << "exp_nrows=" << exp_nRows << " (nVertices= " << nVertices << ", numGhost=" << graph.GetNodeNumGhost() << ")" << std::endl;
228  // std::cout << "nonUniqueMap=" << nonUniqueMap->getNodeNumElements() << std::endl;
229 
230  RCP<Xpetra::Vector<double,LO,GO,NO> > temp_ = Xpetra::VectorFactory<double,LO,GO,NO> ::Build(nonUniqueMap,false); //no need to zero out vector in ctor
231  temp_->putScalar(1.);
232 
233  RCP<Xpetra::Vector<double,LO,GO,NO> > tempOutput_ = Xpetra::VectorFactory<double,LO,GO,NO> ::Build(nonUniqueMap);
234 
235  myWidget.NonUnique2NonUnique(*temp_, *tempOutput_, Xpetra::ADD);
236 
237  std::vector<bool> gidNotShared(exp_nRows);
238  {
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;
243  else
244  gidNotShared[i] = true;
245  }
246  }
247 
248  // Phase 4.
249  double nAggregatesTarget;
250  nAggregatesTarget = ((double) uniqueMap->getGlobalNumElements())* (((double) uniqueMap->getGlobalNumElements())/ ((double) graph.GetGlobalNumEdges()));
251 
252  GO nAggregatesLocal=nAggregates, nAggregatesGlobal; MueLu_sumAll(graph.GetComm(), nAggregatesLocal, nAggregatesGlobal);
253 
254  LO minNAggs; MueLu_minAll(graph.GetComm(), nAggregates, minNAggs);
255  LO maxNAggs; MueLu_maxAll(graph.GetComm(), nAggregates, maxNAggs);
256 
257  //
258  // Only do this phase if things look really bad. THIS
259  // CODE IS PRETTY EXPERIMENTAL
260  //
261 #define MUELU_PHASE4BUCKETS 6
262  if ((nAggregatesGlobal < graph.GetComm()->getSize()) &&
263  (2.5*nAggregatesGlobal < nAggregatesTarget) &&
264  (minNAggs ==0) && (maxNAggs <= 1)) {
265 
266  // Modify seed of the random algorithm used by temp_->randomize()
267  {
268  typedef Teuchos::ScalarTraits<double> scalarTrait; // temp_ is of type double.
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()));
273  }
274 
275  temp_->randomize();
276 
277  ArrayRCP<double> temp = temp_->getDataNonConst(0);
278 
279  // build a list of candidate root nodes (vertices not adjacent
280  // to aggregated vertices)
281 
282  my_size_t nCandidates = 0;
283  global_size_t nCandidatesGlobal;
284 
285  ArrayRCP<LO> candidates = Teuchos::arcp<LO>(nVertices+1);
286 
287  double priorThreshold = 0.;
288  for (int kkk = 0; kkk < MUELU_PHASE4BUCKETS; kkk++) {
289 
290  {
291  ArrayRCP<const LO> vertex2AggId = aggregates.GetVertex2AggId()->getData(0);
292  ArrayView<const LO> vertex2AggIdView = vertex2AggId();
293  RootCandidates(nVertices, vertex2AggIdView, graph, candidates, nCandidates, nCandidatesGlobal);
294  // views on distributed vectors are freed here.
295  }
296 
297  double nTargetNewGuys = nAggregatesTarget - nAggregatesGlobal;
298  double threshold = priorThreshold + (1. - priorThreshold)*nTargetNewGuys/(nCandidatesGlobal + .001);
299 
300  threshold = (threshold*(kkk+1.))/((double) MUELU_PHASE4BUCKETS);
301  priorThreshold = threshold;
302 
303  {
304  ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
305  ArrayRCP<double> weights = distWeights->getDataNonConst(0);
306 
307  for (int k = 0; k < nCandidates; k++ ) {
308  int i = candidates[k];
309  if ((vertex2AggId[i] == MUELU_UNAGGREGATED) && (fabs(temp[i]) < threshold)) {
310  // Note: priorThreshold <= fabs(temp[i]) <= 1
311 
312  // neighOfINode is the neighbor node list of node 'iNode'.
313  ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
314 
315  if (neighOfINode.size() > minNodesPerAggregate) { //TODO: check if this test is exactly was we want to do
316  int count = 0;
317  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
318  int Adjacent = *it;
319  // This might not be true if someone close to i
320  // is chosen as a root via fabs(temp[]) < Threshold
321  if (vertex2AggId[Adjacent] == MUELU_UNAGGREGATED){
322  count++;
323  vertex2AggId[Adjacent] = nAggregates;
324  weights[Adjacent] = 1.;
325  }
326  }
327  if (count >= minNodesPerAggregate) {
328  vertex2AggId[i] = nAggregates++;
329  weights[i] = 2.;
330  aggregates.SetIsRoot(i);
331  }
332  else { // undo things
333  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
334  int Adjacent = *it;
335  if (vertex2AggId[Adjacent] == nAggregates){
336  vertex2AggId[Adjacent] = MUELU_UNAGGREGATED;
337  weights[Adjacent] = 0.;
338  }
339  }
340  }
341  }
342  }
343  }
344  // views on distributed vectors are freed here.
345  }
346  //TODO JJH We want to skip this call
347  myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
348  // All tentatively assigned vertices are now definitive
349  nAggregatesLocal=nAggregates;
350  MueLu_sumAll(graph.GetComm(), nAggregatesLocal, nAggregatesGlobal);
351 
352  // check that there are no aggregates sizes below minNodesPerAggregate
353 
354  aggregates.SetNumAggregates(nAggregates);
355 
356  RemoveSmallAggs(aggregates, minNodesPerAggregate, distWeights, myWidget);
357 
358  nAggregates = aggregates.GetNumAggregates();
359  } // one possibility
360  }
361 
362  // Initialize things for Phase 5. This includes building the transpose
363  // of the matrix ONLY for transposed rows that correspond to unaggregted
364  // ghost vertices. Further, the transpose is only a local transpose.
365  // Nonzero edges which exist on other processors are not represented.
366 
367 
368  int observedNAgg=-1; //number of aggregates that contain vertices on this process
369 
370  {
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];
376  observedNAgg++;
377  }
378 
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);
382 
383  for (int i = 0; i < exp_nRows; i++) Mark[i] = MUELU_DISTONE_VERTEX_WEIGHT;
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;
386 
387  // Grab the transpose matrix graph for unaggregated ghost vertices.
388  // a) count the number of nonzeros per row in the transpose
389  std::vector<int> RowPtr(exp_nRows+1-nVertices);
390  //{
391  ArrayRCP<const LO> vertex2AggIdCst = aggregates.GetVertex2AggId()->getData(0);
392 
393  for (int i = nVertices; i < exp_nRows; i++) RowPtr[i-nVertices] = 0;
394  for (int i = 0; i < nVertices; i++) {
395 
396  // neighOfINode is the neighbor node list of node 'iNode'.
397  ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
398 
399  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
400  int j = *it;
401  if ( (j >= nVertices) && (vertex2AggIdCst[j] == MUELU_UNAGGREGATED)){
402  RowPtr[j-nVertices]++;
403  }
404  }
405  }
406 
407  // b) Convert RowPtr[i] to point to 1st first nnz spot in row i.
408 
409  int iSum = 0, iTemp;
410  for (int i = nVertices; i < exp_nRows; i++) {
411  iTemp = RowPtr[i-nVertices];
412  RowPtr[i-nVertices] = iSum;
413  iSum += iTemp;
414  }
415  RowPtr[exp_nRows-nVertices] = iSum;
416  std::vector<LO> cols(iSum+1);
417 
418  // c) Traverse matrix and insert entries in proper location.
419  for (int i = 0; i < nVertices; i++) {
420 
421  // neighOfINode is the neighbor node list of node 'iNode'.
422  ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
423 
424  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
425  int j = *it;
426  if ( (j >= nVertices) && (vertex2AggIdCst[j] == MUELU_UNAGGREGATED)){
427  cols[RowPtr[j-nVertices]++] = i;
428  }
429  }
430  }
431 
432  // d) RowPtr[i] points to beginning of row i+1 so shift by one location.
433  for (int i = exp_nRows; i > nVertices; i--)
434  RowPtr[i-nVertices] = RowPtr[i-1-nVertices];
435  RowPtr[0] = 0;
436 
437  // views on distributed vectors are freed here.
438  vertex2AggIdCst = Teuchos::null;
439  //}
440 
441  int bestScoreCutoff;
442  int thresholds[10] = {300,200,100,50,25,13,7,4,2,0};
443 
444  // Stick unaggregated vertices into existing aggregates as described above.
445 
446  {
447  int ncalls=0;
448 
449  for (int kk = 0; kk < 10; kk += 2) {
450  bestScoreCutoff = thresholds[kk];
451 
452  ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
453  ArrayRCP<const LO> procWinner = aggregates.GetProcWinner()->getData(0);
454  ArrayRCP<double> weights = distWeights->getDataNonConst(0);
455 
456  for (int i = 0; i < exp_nRows; i++) {
457 
458  if (vertex2AggId[i] == MUELU_UNAGGREGATED) {
459 
460  // neighOfINode is the neighbor node list of node 'iNode'.
461  ArrayView<const LO> neighOfINode;
462 
463  // Grab neighboring vertices which is either in graph for local ids
464  // or sits in transposed fragment just constructed above for ghosts.
465  if (i < nVertices) {
466  neighOfINode = graph.getNeighborVertices(i);
467  }
468  else {
469  LO *rowi_col = NULL, rowi_N;
470  rowi_col = &(cols[RowPtr[i-nVertices]]);
471  rowi_N = RowPtr[i+1-nVertices] - RowPtr[i-nVertices];
472 
473  neighOfINode = ArrayView<const LO>(rowi_col, rowi_N);
474  }
475  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
476  int Adjacent = *it;
477  int AdjacentAgg = vertex2AggId[Adjacent];
478 
479  //Adjacent is aggregated and either I own the aggregate
480  // or I could own the aggregate after arbitration.
481  if ((AdjacentAgg != MUELU_UNAGGREGATED) &&
482  ((procWinner[Adjacent] == myPid) ||
483  (procWinner[Adjacent] == MUELU_UNASSIGNED))){
484  SumOfMarks[AdjacentAgg] += Mark[Adjacent];
485  }
486  }
487  int best_score = MUELU_NOSCORE;
488  int best_agg = -1;
489  int BestMark = -1;
490  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).
491 
492  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
493  int Adjacent = *it;
494  int AdjacentAgg = vertex2AggId[Adjacent];
495  //Adjacent is unaggregated, has some value and no
496  //other processor has definitively claimed him
497  if ((AdjacentAgg != MUELU_UNAGGREGATED) &&
498  (SumOfMarks[AdjacentAgg] != 0) &&
499  ((procWinner[Adjacent] == myPid) ||
500  (procWinner[Adjacent] == MUELU_UNASSIGNED ))) {
501 
502  // first figure out the penalty associated with
503  // AdjacentAgg having already been incremented
504  // during this phase, then compute score.
505 
506  double penalty = (double) (INCR_SCALING*agg_incremented[AdjacentAgg]);
507  if (penalty > MUELU_PENALTYFACTOR*((double)SumOfMarks[AdjacentAgg]))
508  penalty = MUELU_PENALTYFACTOR*((double)SumOfMarks[AdjacentAgg]);
509  int score = SumOfMarks[AdjacentAgg]- ((int) floor(penalty));
510 
511  if (score > best_score) {
512  best_agg = AdjacentAgg;
513  best_score = score;
514  BestMark = Mark[Adjacent];
515  cannotLoseAllFriends = false;
516 
517  // This address issue mentioned above by checking whether
518  // Adjacent could be lost in arbitration. weight==0 means that
519  // Adjacent was not set during this loop of Phase 5 (and so it
520  // has already undergone arbitration). GidNotShared == true
521  // obviously implies that Adjacent cannot be lost to arbitration
522  if ((weights[Adjacent]== 0.) || (gidNotShared[Adjacent] == true))
523  cannotLoseAllFriends = true;
524  }
525  // Another vertex within current best aggregate found.
526  // We should have (best_score == score). We need to see
527  // if we can improve BestMark and cannotLoseAllFriends.
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];
532  }
533  }
534  }
535  // Clean up
536  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
537  int Adjacent = *it;
538  int AdjacentAgg = vertex2AggId[Adjacent];
539  if (AdjacentAgg >= 0) SumOfMarks[AdjacentAgg] = 0;
540  }
541  // Tentatively assign vertex to best_agg.
542  if ( (best_score >= bestScoreCutoff) && (cannotLoseAllFriends)) {
543 
544  TEUCHOS_TEST_FOR_EXCEPTION(best_agg == -1 || BestMark == -1, MueLu::Exceptions::RuntimeError, "MueLu::CoupledAggregationFactory internal error"); // should never happen
545 
546  vertex2AggId[i] = best_agg;
547  weights[i] = best_score;
548  agg_incremented[best_agg]++;
549  Mark[i] = (int) ceil( ((double) BestMark)/2.);
550  }
551  }
552 
553  // views on distributed vectors are freed here.
554  }
555 
556  vertex2AggId = Teuchos::null;
557  procWinner = Teuchos::null;
558  weights = Teuchos::null;
559 
560  ++ncalls;
561  //TODO JJH We want to skip this call
562  myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
563  // All tentatively assigned vertices are now definitive
564  }
565 
566  // if (graph.GetComm()->getRank()==0)
567  // std::cout << "#calls to Arb&Comm=" << ncalls << std::endl;
568  }
569 
570  // Phase 6: Aggregate remain unaggregated vertices and try at all costs
571  // to avoid small aggregates.
572  // One case where we can find ourselves in this situation
573  // is if all vertices vk adjacent to v have already been
574  // put in other processor's aggregates and v does not have
575  // a direct connection to a local vertex in any of these
576  // aggregates.
577 
578  int Nleftover = 0, Nsingle = 0;
579  {
580 
581  ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
582  ArrayRCP<double> weights = distWeights->getDataNonConst(0);
583  ArrayRCP<const LO> procWinner = aggregates.GetProcWinner()->getData(0);
584 
585  int count = 0;
586  for (my_size_t i = 0; i < nVertices; i++) {
587  if (vertex2AggId[i] == MUELU_UNAGGREGATED) {
588  Nleftover++;
589 
590  // neighOfINode is the neighbor node list of node 'iNode'.
591  ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
592 
593  // We don't want too small of an aggregate. So lets see if there is an
594  // unaggregated neighbor that we can also put with this vertex
595 
596  vertex2AggId[i] = nAggregates;
597  weights[i] = 1.;
598  if (count == 0) aggregates.SetIsRoot(i);
599  count++;
600  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
601  int j = *it;
602  if ((j != i)&&(vertex2AggId[j] == MUELU_UNAGGREGATED)&&
603  (j < nVertices)) {
604  vertex2AggId[j] = nAggregates;
605  weights[j] = 1.;
606  count++;
607  }
608  }
609  if ( count >= minNodesPerAggregate) {
610  nAggregates++;
611  count = 0;
612  }
613  }
614  }
615 
616  // We have something which is under minNodesPerAggregate when
617  if (count != 0) {
618 #ifdef FIXME
619  // Can stick small aggregate with 0th aggregate?
620  if (nAggregates > 0) {
621  for (my_size_t i = 0; i < nVertices; i++) {
622  if ((vertex2AggId[i] == nAggregates) && (procWinner[i] == myPid)) {
623  vertex2AggId[i] = 0;
624  aggregates.SetIsRoot(i,false);
625  }
626  }
627  }
628  else {
629  Nsingle++;
630  nAggregates++;
631  }
632 #else
633  // Can stick small aggregate with 0th aggregate?
634  if (nAggregates > 0) {
635  for (my_size_t i = 0; i < nVertices; i++) {
636  // TW: This is not a real fix. This may produce ugly bad aggregates!
637  // I removed the procWinner[i] == myPid check. it makes no sense to me since
638  // it leaves vertex2AggId[i] == nAggregates -> crash in ComputeAggregateSizes().
639  // Maybe it's better to add the leftovers to the last generated agg on the current proc.
640  // The best solution would be to add them to the "next"/nearest aggregate, that may be
641  // on an other processor
642  if (vertex2AggId[i] == nAggregates) {
643  vertex2AggId[i] = nAggregates-1; //0;
644  aggregates.SetIsRoot(i,false);
645  }
646  }
647  }
648  else {
649  Nsingle++;
650  nAggregates++;
651  }
652 #endif
653  }
654 
655  // views on distributed vectors are freed here.
656  }
657 
658  //TODO JJH We want to skip this call
659  myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, false);
660 
661  if (IsPrint(Statistics1)) {
662  GO total_Nsingle=0; MueLu_sumAll(graph.GetComm(), (GO)Nsingle, total_Nsingle);
663  GO total_Nleftover=0; MueLu_sumAll(graph.GetComm(), (GO)Nleftover, total_Nleftover);
664  // GO total_aggs; MueLu_sumAll(graph.GetComm(), (GO)nAggregates, total_aggs);
665  // GetOStream(Statistics1) << "Phase 6 - total aggregates = " << total_aggs << std::endl;
666  GetOStream(Statistics1) << "Phase 6 - leftovers = " << total_Nleftover << " and singletons = " << total_Nsingle << std::endl;
667  }
668 
669  aggregates.SetNumAggregates(nAggregates);
670 
671  } //AggregateLeftovers
672 
673  template <class LocalOrdinal, class GlobalOrdinal, class Node>
675  ArrayView<const LO> & vertex2AggId, GraphBase const &graph,
676  ArrayRCP<LO> &candidates, my_size_t &nCandidates, global_size_t &nCandidatesGlobal) const
677  {
678  nCandidates = 0;
679 
680  for (my_size_t i = 0; i < nVertices; i++ ) {
681  if (vertex2AggId[i] == MUELU_UNAGGREGATED) {
682  bool noAggdNeighbors = true;
683 
684  // neighOfINode is the neighbor node list of node 'iNode'.
685  ArrayView<const LO> neighOfINode = graph.getNeighborVertices(i);
686 
687  for (typename ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it) {
688  int adjacent = *it;
689  if (vertex2AggId[adjacent] != MUELU_UNAGGREGATED)
690  noAggdNeighbors = false;
691  }
692  if (noAggdNeighbors == true) candidates[nCandidates++] = i;
693  }
694  }
695 
696  MueLu_sumAll(graph.GetComm(), (GO)nCandidates, nCandidatesGlobal);
697 
698  } //RootCandidates
699 
700  template <class LocalOrdinal, class GlobalOrdinal, class Node>
702  RCP<Xpetra::Vector<double,LO,GO,NO> > & distWeights, const MueLu::CoupledAggregationCommHelper<LO,GO,NO> & myWidget) const {
703  int myPid = aggregates.GetMap()->getComm()->getRank();
704 
705  LO nAggregates = aggregates.GetNumAggregates();
706 
707  ArrayRCP<LO> procWinner = aggregates.GetProcWinner()->getDataNonConst(0);
708  ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
709  LO size = procWinner.size();
710 
711  //ArrayRCP<int> AggInfo = Teuchos::arcp<int>(nAggregates+1);
712  //aggregates.ComputeAggSizes(AggInfo);
713  ArrayRCP<LO> AggInfo = aggregates.ComputeAggregateSizes();
714 
715  ArrayRCP<double> weights = distWeights->getDataNonConst(0);
716 
717  // Make a list of all aggregates indicating New AggId
718  // Use AggInfo array for this.
719 
720  LO NewNAggs = 0;
721  for (LO i = 0; i < nAggregates; i++) {
722  if ( AggInfo[i] < min_size) {
723  AggInfo[i] = MUELU_UNAGGREGATED;
724  }
725  else AggInfo[i] = NewNAggs++;
726  }
727 
728  for (LO k = 0; k < size; k++ ) {
729  if (procWinner[k] == myPid) {
730  if (vertex2AggId[k] != MUELU_UNAGGREGATED) {
731  vertex2AggId[k] = AggInfo[vertex2AggId[k]];
732  weights[k] = 1.;
733  }
734  if (vertex2AggId[k] == MUELU_UNAGGREGATED)
735  aggregates.SetIsRoot(k,false);
736  }
737  }
738  nAggregates = NewNAggs;
739 
740  //TODO JJH We want to skip this call
741  myWidget.ArbitrateAndCommunicate(*distWeights, aggregates, true);
742  // All tentatively assigned vertices are now definitive
743 
744  // procWinner is not set correctly for aggregates which have
745  // been eliminated
746  for (LO i = 0; i < size; i++) {
747  if (vertex2AggId[i] == MUELU_UNAGGREGATED)
748  procWinner[i] = MUELU_UNASSIGNED;
749  }
750  aggregates.SetNumAggregates(nAggregates);
751 
752  return 0; //TODO
753  } //RemoveSmallAggs
754 
755 } //namespace MueLu
756 
757 #endif // MUELU_LEFTOVERAGGREGATIONALGORITHM_DEF_HPP
#define MUELU_UNASSIGNED
#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.
Print more statistics.
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
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.
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
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.
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 &#39;v&#39;.
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.