MueLu  Version of the Day
MueLu_LocalAggregationAlgorithm_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_LOCALAGGREGATIONALGORITHM_DEF_HPP
47 #define MUELU_LOCALAGGREGATIONALGORITHM_DEF_HPP
48 
49 #include <Teuchos_Comm.hpp>
50 #include <Teuchos_CommHelpers.hpp>
51 
52 #include <Xpetra_Vector.hpp>
53 
55 
56 #include "MueLu_Aggregates.hpp"
57 #include "MueLu_Exceptions.hpp"
58 #include "MueLu_GraphBase.hpp"
59 #include "MueLu_LinkedList.hpp"
60 #include "MueLu_Monitor.hpp"
61 #include "MueLu_Utilities.hpp"
62 
63 namespace MueLu {
64 
65  template <class LocalOrdinal, class GlobalOrdinal, class Node>
67  : ordering_("natural"), minNodesPerAggregate_(1), maxNeighAlreadySelected_(0)
68  { }
69 
70  template <class LocalOrdinal, class GlobalOrdinal, class Node>
72  Monitor m(*this, "Coarsen Uncoupled");
73 
74  GetOStream(Runtime1) << "Ordering: " << ordering_ << std::endl;
75  GetOStream(Runtime1) << "Min nodes per aggregate: " << minNodesPerAggregate_ << std::endl;
76  GetOStream(Runtime1) << "Max nbrs already selected: " << maxNeighAlreadySelected_ << std::endl;
77 
78  /* Create Aggregation object */
79  my_size_t nAggregates = 0;
80 
81  /* ============================================================= */
82  /* aggStat indicates whether this node has been aggreated, and */
83  /* vertex2AggId stores the aggregate number where this node has */
84  /* been aggregated into. */
85  /* ============================================================= */
86 
87  Teuchos::ArrayRCP<CANodeState> aggStat;
88  const my_size_t nRows = graph.GetNodeNumVertices();
89  if (nRows > 0) aggStat = Teuchos::arcp<CANodeState>(nRows);
90  for ( my_size_t i = 0; i < nRows; ++i ) aggStat[i] = CA_READY;
91 
92  /* ============================================================= */
93  /* Phase 1 : */
94  /* for all nodes, form a new aggregate with its neighbors */
95  /* if the number of its neighbors having been aggregated does */
96  /* not exceed a given threshold */
97  /* (GetMaxNeighAlreadySelected() = 0 ===> Vanek's scheme) */
98  /* ============================================================= */
99 
100  /* some general variable declarations */
101  Teuchos::ArrayRCP<LO> randomVector;
102  RCP<MueLu::LinkedList> nodeList; /* list storing the next node to pick as a root point for ordering_ == "graph" */
103  MueLu_SuperNode *aggHead=NULL, *aggCurrent=NULL, *supernode=NULL;
104 
105 
106  if ( ordering_ == "random" ) /* random ordering */
107  {
108  //TODO: could be stored in a class that respect interface of LinkedList
109 
110  randomVector = Teuchos::arcp<LO>(nRows); //size_t or int ?-> to be propagated
111  for (my_size_t i = 0; i < nRows; ++i) randomVector[i] = i;
112  RandomReorder(randomVector);
113  }
114  else if ( ordering_ == "graph" ) /* graph ordering */
115  {
116  nodeList = rcp(new MueLu::LinkedList());
117  nodeList->Add(0);
118  }
119 
120  /* main loop */
121  {
122  LO iNode = 0;
123  LO iNode2 = 0;
124 
125  Teuchos::ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0); // output only: contents ignored
126 
127  while (iNode2 < nRows)
128  {
129 
130  /*------------------------------------------------------ */
131  /* pick the next node to aggregate */
132  /*------------------------------------------------------ */
133 
134  if ( ordering_ == "natural" ) iNode = iNode2++;
135  else if ( ordering_ == "random" ) iNode = randomVector[iNode2++];
136  else if ( ordering_ == "graph" )
137  {
138  if ( nodeList->IsEmpty() )
139  {
140  for ( int jNode = 0; jNode < nRows; ++jNode )
141  {
142  if ( aggStat[jNode] == CA_READY )
143  {
144  nodeList->Add(jNode); //TODO optim: not necessary to create a node. Can just set iNode value and skip the end
145  break;
146  }
147  }
148  }
149  if ( nodeList->IsEmpty() ) break; /* end of the while loop */ //TODO: coding style :(
150 
151  iNode = nodeList->Pop();
152  }
153  else {
154  throw(Exceptions::RuntimeError("CoarsenUncoupled: bad aggregation ordering option"));
155  }
156 
157  /*------------------------------------------------------ */
158  /* consider further only if the node is in CA_READY mode */
159  /*------------------------------------------------------ */
160 
161  if ( aggStat[iNode] == CA_READY )
162  {
163  // neighOfINode is the neighbor node list of node 'iNode'.
164  Teuchos::ArrayView<const LO> neighOfINode = graph.getNeighborVertices(iNode);
165  typename Teuchos::ArrayView<const LO>::size_type length = neighOfINode.size();
166 
167  supernode = new MueLu_SuperNode;
168  try {
169  supernode->list = Teuchos::arcp<int>(length+1);
170  } catch (std::bad_alloc&) {
171  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::LocalAggregationAlgorithm::CoarsenUncoupled(): Error: couldn't allocate memory for supernode! length=" + Teuchos::toString(length));
172  }
173 
174  supernode->maxLength = length;
175  supernode->length = 1;
176  supernode->list[0] = iNode;
177 
178  int selectFlag = 1;
179  {
180  /*--------------------------------------------------- */
181  /* count the no. of neighbors having been aggregated */
182  /*--------------------------------------------------- */
183 
184  int count = 0;
185  for (typename Teuchos::ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it)
186  {
187  int index = *it;
188  if ( index < nRows )
189  {
190  if ( aggStat[index] == CA_READY ||
191  aggStat[index] == CA_NOTSEL )
192  supernode->list[supernode->length++] = index;
193  else count++;
194 
195  }
196  }
197 
198  /*--------------------------------------------------- */
199  /* if there are too many neighbors aggregated or the */
200  /* number of nodes in the new aggregate is too few, */
201  /* don't do this one */
202  /*--------------------------------------------------- */
203 
204  if ( count > GetMaxNeighAlreadySelected() ) selectFlag = 0;
205  }
206 
207  // Note: the supernode length is actually 1 more than the
208  // number of nodes in the candidate aggregate. The
209  // root is counted twice. I'm not sure if this is
210  // a bug or a feature ... so I'll leave it and change
211  // < to <= in the if just below.
212 
213  if (selectFlag != 1 ||
214  supernode->length <= GetMinNodesPerAggregate())
215  {
216  aggStat[iNode] = CA_NOTSEL;
217  delete supernode;
218  if ( ordering_ == "graph" ) /* if graph ordering */
219  {
220  for (typename Teuchos::ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it)
221  {
222  int index = *it;
223  if ( index < nRows && aggStat[index] == CA_READY )
224  {
225  nodeList->Add(index);
226  }
227  }
228  }
229  }
230  else
231  {
232  aggregates.SetIsRoot(iNode);
233  for ( int j = 0; j < supernode->length; ++j )
234  {
235  int jNode = supernode->list[j];
236  aggStat[jNode] = CA_SELECTED;
237  vertex2AggId[jNode] = nAggregates;
238  if ( ordering_ == "graph" ) /* if graph ordering */
239  {
240 
241  Teuchos::ArrayView<const LO> neighOfJNode = graph.getNeighborVertices(jNode);
242 
243  for (typename Teuchos::ArrayView<const LO>::const_iterator it = neighOfJNode.begin(); it != neighOfJNode.end(); ++it)
244  {
245  int index = *it;
246  if ( index < nRows && aggStat[index] == CA_READY )
247  {
248  nodeList->Add(index);
249  }
250  }
251  }
252  }
253  supernode->next = NULL;
254  supernode->index = nAggregates;
255  if ( nAggregates == 0 )
256  {
257  aggHead = supernode;
258  aggCurrent = supernode;
259  }
260  else
261  {
262  aggCurrent->next = supernode;
263  aggCurrent = supernode;
264  }
265  nAggregates++;
266  // unused aggCntArray[nAggregates] = supernode->length;
267  }
268  }
269  } // end of 'for'
270 
271  // views on distributed vectors are freed here.
272 
273  } // end of 'main loop'
274 
275  nodeList = Teuchos::null;
276 
277  /* Update aggregate object */
278  aggregates.SetNumAggregates(nAggregates);
279 
280  /* Verbose */
281  {
282  const RCP<const Teuchos::Comm<int> > & comm = graph.GetComm();
283 
284  if (IsPrint(Warnings0)) {
285  GO localReady=0, globalReady;
286 
287  // Compute 'localReady'
288  for ( my_size_t i = 0; i < nRows; ++i )
289  if (aggStat[i] == CA_READY) localReady++;
290 
291  // Compute 'globalReady'
292  MueLu_sumAll(comm, localReady, globalReady);
293 
294  if(globalReady > 0)
295  GetOStream(Warnings0) << globalReady << " CA_READY nodes left" << std::endl;
296  }
297 
298  if (IsPrint(Statistics1)) {
299  // Compute 'localSelected'
300  LO localSelected=0;
301  for ( my_size_t i = 0; i < nRows; ++i )
302  if ( aggStat[i] == CA_SELECTED ) localSelected++;
303 
304  // Compute 'globalSelected'
305  GO globalSelected; MueLu_sumAll(comm, (GO)localSelected, globalSelected);
306 
307  // Compute 'globalNRows'
308  GO globalNRows; MueLu_sumAll(comm, (GO)nRows, globalNRows);
309 
310  GetOStream(Statistics1) << "Nodes aggregated = " << globalSelected << " (" << globalNRows << ")" << std::endl;
311  }
312 
313  if (IsPrint(Statistics1)) {
314  GO nAggregatesGlobal; MueLu_sumAll(comm, (GO)nAggregates, nAggregatesGlobal);
315  GetOStream(Statistics1) << "Total aggregates = " << nAggregatesGlobal << std::endl;
316  }
317 
318  } // verbose
319 
320  /* ------------------------------------------------------------- */
321  /* clean up */
322  /* ------------------------------------------------------------- */
323 
324  aggCurrent = aggHead;
325  while ( aggCurrent != NULL )
326  {
327  supernode = aggCurrent;
328  aggCurrent = aggCurrent->next;
329  delete supernode;
330  }
331 
332  } // CoarsenUncoupled
333 
334  template <class LocalOrdinal, class GlobalOrdinal, class Node>
336  //TODO: replace int
337  int n = list.size();
338  for(int i=0; i<n-1; i++) {
339  std::swap(list[i], list[RandomOrdinal(i,n-1)]);
340  }
341  }
342 
343  template <class LocalOrdinal, class GlobalOrdinal, class Node>
345  return min + Teuchos::as<int>((max-min+1) * (static_cast<double>(std::rand()) / (RAND_MAX + 1.0)));
346  }
347 
348 } // namespace MueLu
349 
350 #endif // MUELU_LOCALAGGREGATIONALGORITHM_DEF_HPP
Important warning messages (one line)
LocalOrdinal RandomOrdinal(LocalOrdinal min, LocalOrdinal max)
#define MueLu_sumAll(rcpComm, in, out)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Container class for aggregation information.
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
Print more statistics.
Namespace for MueLu classes and methods.
void SetIsRoot(LO i, bool value=true)
Set root node information.
void RandomReorder(Teuchos::ArrayRCP< LO > list) const
Utility to take a list of integers and reorder them randomly (by using a local permutation).
virtual const RCP< const Teuchos::Comm< int > > GetComm() const =0
MueLu representation of a graph.
struct MueLu::MueLu_SuperNode_Struct MueLu_SuperNode
Timer to be used in non-factories.
void RandomReorder(Teuchos::Array< LocalOrdinal > &list)
int RandomOrdinal(int min, int max) const
Generate a random number in the range [min, max].
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
void CoarsenUncoupled(GraphBase const &graph, Aggregates &aggregates) const
Local aggregation.
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.