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