Zoltan2
Zoltan2_AlgRCM.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Zoltan2: A package of combinatorial algorithms for scientific computing
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 Karen Devine (kddevin@sandia.gov)
39 // Erik Boman (egboman@sandia.gov)
40 // Siva Rajamanickam (srajama@sandia.gov)
41 //
42 // ***********************************************************************
43 //
44 // @HEADER
45 #ifndef _ZOLTAN2_ALGRCM_HPP_
46 #define _ZOLTAN2_ALGRCM_HPP_
47 
48 #include <Zoltan2_Algorithm.hpp>
49 #include <Zoltan2_GraphModel.hpp>
51 #include <Zoltan2_Sort.hpp>
52 #include <queue>
53 
54 
58 
59 
60 namespace Zoltan2{
61 
62 template <typename Adapter>
63 class AlgRCM : public Algorithm<Adapter>
64 {
65  private:
66 
67  const RCP<GraphModel<Adapter> > model;
68  const RCP<Teuchos::ParameterList> pl;
69  const RCP<const Teuchos::Comm<int> > comm;
70 
71  public:
72 
73  typedef typename Adapter::lno_t lno_t;
74  typedef typename Adapter::gno_t gno_t;
75  typedef typename Adapter::scalar_t scalar_t;
76 
78  const RCP<GraphModel<Adapter> > &model__,
79  const RCP<Teuchos::ParameterList> &pl__,
80  const RCP<const Teuchos::Comm<int> > &comm__
81  ) : model(model__), pl(pl__), comm(comm__)
82  {
83  }
84 
85  int order(const RCP<OrderingSolution<lno_t, gno_t> > &solution)
86  {
87  int ierr= 0;
88 
89  HELLO;
90 
91  // Get local graph.
92  ArrayView<const gno_t> edgeIds;
93  ArrayView<const lno_t> offsets;
94  ArrayView<StridedData<lno_t, scalar_t> > wgts;
95 
96  const size_t nVtx = model->getLocalNumVertices();
97  model->getEdgeList(edgeIds, offsets, wgts);
98  const int numWeightsPerEdge = model->getNumWeightsPerEdge();
99  if (numWeightsPerEdge > 1){
100  throw std::runtime_error("Multiple weights not supported.");
101  }
102 
103 #if 0
104  // Debug
105  cout << "Debug: Local graph from getLocalEdgeList" << endl;
106  cout << "rank " << comm->getRank() << ": nVtx= " << nVtx << endl;
107  cout << "rank " << comm->getRank() << ": edgeIds: " << edgeIds << endl;
108  cout << "rank " << comm->getRank() << ": offsets: " << offsets << endl;
109 #endif
110 
111  // RCM constructs invPerm, not perm
112  const ArrayRCP<lno_t> invPerm = solution->getPermutationRCP(true);
113 
114  // Check if there are actually edges to reorder.
115  // If there are not, then just use the natural ordering.
116  if (offsets[nVtx] == 0) {
117  for (size_t i = 0; i < nVtx; ++i) {
118  invPerm[i] = i;
119  }
120  solution->setHaveInverse(true);
121  return 0;
122  }
123 
124  // Set the label of each vertex to invalid.
125  Tpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
126  for (size_t i = 0; i < nVtx; ++i) {
127  invPerm[i] = INVALID;
128  }
129 
130  // Loop over all connected components.
131  // Do BFS within each component.
132  gno_t root = 0;
133  std::queue<gno_t> Q;
134  size_t count = 0; // CM label, reversed later
135  size_t next = 0; // next unmarked vertex
136  Teuchos::Array<std::pair<gno_t, size_t> > children; // children and their degrees
137 
138  while (count < nVtx) {
139 
140  // Find suitable root vertex for this component.
141  // First find an unmarked vertex, use to find root in next component.
142  while ((next < nVtx) && (static_cast<Tpetra::global_size_t>(invPerm[next]) != INVALID)) next++;
143 
144  // Select root method. Pseudoperipheral usually gives the best
145  // ordering, but the user may choose a faster method.
146  std::string root_method = pl->get("root_method", "pseudoperipheral");
147  if (root_method == std::string("first"))
148  root = next;
149  else if (root_method == std::string("smallest_degree"))
150  root = findSmallestDegree(next, nVtx, edgeIds, offsets);
151  else if (root_method == std::string("pseudoperipheral"))
152  root = findPseudoPeripheral(next, nVtx, edgeIds, offsets);
153  else {
154  // This should never happen if pl was validated.
155  throw std::runtime_error("invalid root_method");
156  }
157 
158  // Label connected component starting at root
159  Q.push(root);
160  //cout << "Debug: invPerm[" << root << "] = " << count << endl;
161  invPerm[root] = count++;
162 
163  while (Q.size()){
164  // Get a vertex from the queue
165  gno_t v = Q.front();
166  Q.pop();
167  //cout << "Debug: v= " << v << ", offsets[v] = " << offsets[v] << endl;
168 
169  // Add unmarked children to list of pairs, to be added to queue.
170  children.resize(0);
171  for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
172  gno_t child = edgeIds[ptr];
173  if (static_cast<Tpetra::global_size_t>(invPerm[child]) == INVALID){
174  // Not visited yet; add child to list of pairs.
175  std::pair<gno_t,size_t> newchild;
176  newchild.first = child;
177  newchild.second = offsets[child+1] - offsets[child];
178  children.push_back(newchild);
179  }
180  }
181  // Sort children by increasing degree
182  // TODO: If edge weights, sort children by decreasing weight,
184  zort.sort(children);
185 
186  typename Teuchos::Array<std::pair<gno_t,size_t> >::iterator it = children.begin ();
187  for ( ; it != children.end(); ++it){
188  // Push children on the queue in sorted order.
189  gno_t child = it->first;
190  invPerm[child] = count++; // Label as we push on Q
191  Q.push(child);
192  //cout << "Debug: invPerm[" << child << "] = " << count << endl;
193  }
194  }
195  }
196 
197  // Reverse labels for RCM
198  bool reverse = true; // TODO: Make parameter
199  if (reverse) {
200  lno_t temp;
201  for (size_t i=0; i < nVtx/2; ++i) {
202  // Swap (invPerm[i], invPerm[nVtx-i])
203  temp = invPerm[i];
204  invPerm[i] = invPerm[nVtx-1-i];
205  invPerm[nVtx-1-i] = temp;
206  }
207  }
208 
209  solution->setHaveInverse(true);
210  return ierr;
211  }
212 
213  private:
214  // Find a smallest degree vertex in component containing v
215  gno_t findSmallestDegree(
216  gno_t v,
217  lno_t nVtx,
218  ArrayView<const gno_t> edgeIds,
219  ArrayView<const lno_t> offsets)
220  {
221  std::queue<gno_t> Q;
222  Teuchos::Array<bool> mark(nVtx);
223 
224  // Do BFS and compute smallest degree as we go
225  lno_t smallestDegree = nVtx;
226  gno_t smallestVertex = 0;
227 
228  // Clear mark array - nothing marked yet
229  for (int i=0; i<nVtx; i++)
230  mark[i] = false;
231 
232  // Start from v
233  Q.push(v);
234  while (Q.size()){
235  // Get first vertex from the queue
236  v = Q.front();
237  Q.pop();
238  // Check degree of v
239  lno_t deg = offsets[v+1] - offsets[v];
240  if (deg < smallestDegree){
241  smallestDegree = deg;
242  smallestVertex = v;
243  }
244  // Add unmarked children to queue
245  for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
246  gno_t child = edgeIds[ptr];
247  if (!mark[child]){
248  mark[child] = true;
249  Q.push(child);
250  }
251  }
252  }
253  return smallestVertex;
254  }
255 
256  // Find a pseudoperipheral vertex in component containing v
257  gno_t findPseudoPeripheral(
258  gno_t v,
259  lno_t nVtx,
260  ArrayView<const gno_t> edgeIds,
261  ArrayView<const lno_t> offsets)
262  {
263  std::queue<gno_t> Q;
264  Teuchos::Array<bool> mark(nVtx);
265 
266  // Do BFS a couple times, pick vertex last visited (furthest away)
267  const int numBFS = 2;
268  for (int bfs=0; bfs<numBFS; bfs++){
269  // Clear mark array - nothing marked yet
270  for (int i=0; i<nVtx; i++)
271  mark[i] = false;
272  // Start from v
273  Q.push(v);
274  while (Q.size()){
275  // Get first vertex from the queue
276  v = Q.front();
277  Q.pop();
278  // Add unmarked children to queue
279  for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
280  gno_t child = edgeIds[ptr];
281  if (!mark[child]){
282  mark[child] = true;
283  Q.push(child);
284  }
285  }
286  }
287  }
288  return v;
289  }
290 
291 };
292 }
293 #endif
#define HELLO
#define INVALID(STR)
size_t global_size_t
Defines the OrderingSolution class.
AlgRCM(const RCP< GraphModel< Adapter > > &model__, const RCP< Teuchos::ParameterList > &pl__, const RCP< const Teuchos::Comm< int > > &comm__)
int order(const RCP< OrderingSolution< lno_t, gno_t > > &solution)
Ordering method.
Algorithm defines the base class for all algorithms.
Sort vector of pairs (key, value) by value.
Adapter::scalar_t scalar_t
GraphModel defines the interface required for graph models.
Defines the GraphModel interface.
The class containing ordering solutions.
Adapter::lno_t lno_t
void sort(Teuchos::Array< std::pair< key_t, value_t > > &listofPairs, bool inc=true)
Adapter::gno_t gno_t