MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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)
#define MueLu_sumAll(rcpComm, in, out)
Container class for aggregation information.
void CoarsenUncoupled(GraphBase const &graph, Aggregates &aggregates) const
Local aggregation.
GlobalOrdinal GO
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)
const_pointer const_iterator
Print more statistics.
size_type size() const
LocalOrdinal LO
size_type size() const
void RandomReorder(Teuchos::ArrayRCP< LO > list) const
Utility to take a list of integers and reorder them randomly (by using a local permutation).
void SetIsRoot(LO i, bool value=true)
Set root node information.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
virtual const RCP< const Teuchos::Comm< int > > GetComm() const =0
MueLu representation of a graph.
void Add(int iNode)
struct MueLu::MueLu_SuperNode_Struct MueLu_SuperNode
Timer to be used in non-factories.
iterator end() const
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
virtual Teuchos::ArrayView< const LocalOrdinal > getNeighborVertices(LocalOrdinal v) const =0
Return the list of vertices adjacent to the vertex &#39;v&#39;.
std::string toString(const T &t)
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.