Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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<const typename Adapter::base_adapter_t> adapter;
68  const RCP<Teuchos::ParameterList> pl;
69  const RCP<const Teuchos::Comm<int> > comm;
70  RCP<const Environment> env;
71  modelFlag_t graphFlags;
72 
73  public:
74 
75  typedef typename Adapter::lno_t lno_t;
76  typedef typename Adapter::gno_t gno_t;
77  typedef typename Adapter::offset_t offset_t;
78  typedef typename Adapter::scalar_t scalar_t;
79 
81  const RCP<const typename Adapter::base_adapter_t> &adapter__,
82  const RCP<Teuchos::ParameterList> &pl__,
83  const RCP<const Teuchos::Comm<int> > &comm__,
84  RCP<const Environment> &env__,
85  const modelFlag_t &graphFlags__
86  ) : adapter(adapter__), pl(pl__), comm(comm__), env(env__), graphFlags(graphFlags__)
87  {
88  }
89 
90  int globalOrder(const RCP<GlobalOrderingSolution<gno_t> > &/* solution */)
91  {
92  throw std::logic_error("AlgRCM does not yet support global ordering.");
93  }
94 
95  int localOrder(const RCP<LocalOrderingSolution<lno_t> > &solution)
96  {
97  int ierr= 0;
98 
99  HELLO;
100 
101  // Get local graph.
102  ArrayView<const gno_t> edgeIds;
103  ArrayView<const offset_t> offsets;
104  ArrayView<StridedData<lno_t, scalar_t> > wgts;
105 
106  const auto model = rcp(new GraphModel<Adapter>(adapter, env, comm, graphFlags));
107  const size_t nVtx = model->getLocalNumVertices();
108  model->getEdgeList(edgeIds, offsets, wgts);
109  const int numWeightsPerEdge = model->getNumWeightsPerEdge();
110  if (numWeightsPerEdge > 1){
111  throw std::runtime_error("Multiple weights not supported.");
112  }
113 
114 #if 0
115  // Debug
116  cout << "Debug: Local graph from getLocalEdgeList" << endl;
117  cout << "rank " << comm->getRank() << ": nVtx= " << nVtx << endl;
118  cout << "rank " << comm->getRank() << ": edgeIds: " << edgeIds << endl;
119  cout << "rank " << comm->getRank() << ": offsets: " << offsets << endl;
120 #endif
121 
122  // RCM constructs invPerm, not perm
123  const ArrayRCP<lno_t> invPerm = solution->getPermutationRCP(true);
124  const ArrayRCP<lno_t> tmpPerm(invPerm.size()); //temporary array used in reversing order
125 
126  // Check if there are actually edges to reorder.
127  // If there are not, then just use the natural ordering.
128  if (offsets[nVtx] == 0) {
129  for (size_t i = 0; i < nVtx; ++i) {
130  invPerm[i] = i;
131  }
132  solution->setHaveInverse(true);
133  return 0;
134  }
135 
136  // Set the label of each vertex to invalid.
137  Tpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
138  for (size_t i = 0; i < nVtx; ++i) {
139  invPerm[i] = INVALID;
140  }
141 
142  // Loop over all connected components.
143  // Do BFS within each component.
144  gno_t root = 0;
145  std::queue<gno_t> Q;
146  size_t count = 0; // CM label, reversed later
147  size_t next = 0; // next unmarked vertex
148  Teuchos::Array<std::pair<gno_t, offset_t> > children; // children and their degrees
149 
150  while (count < nVtx) {
151 
152  // Find suitable root vertex for this component.
153  // First find an unmarked vertex, use to find root in next component.
154  while ((next < nVtx) && (static_cast<Tpetra::global_size_t>(invPerm[next]) != INVALID)) next++;
155 
156  // Select root method. Pseudoperipheral usually gives the best
157  // ordering, but the user may choose a faster method.
158  std::string root_method = pl->get("root_method", "pseudoperipheral");
159  if (root_method == std::string("first"))
160  root = next;
161  else if (root_method == std::string("smallest_degree"))
162  root = findSmallestDegree(next, nVtx, edgeIds, offsets);
163  else if (root_method == std::string("pseudoperipheral"))
164  root = findPseudoPeripheral(next, nVtx, edgeIds, offsets);
165  else {
166  // This should never happen if pl was validated.
167  throw std::runtime_error("invalid root_method");
168  }
169 
170  // Label connected component starting at root
171  Q.push(root);
172  //cout << "Debug: invPerm[" << root << "] = " << count << endl;
173  invPerm[root] = count++;
174  tmpPerm[invPerm[root]] = root;
175 
176  while (Q.size()){
177  // Get a vertex from the queue
178  gno_t v = Q.front();
179  Q.pop();
180  //cout << "Debug: v= " << v << ", offsets[v] = " << offsets[v] << endl;
181 
182  // Add unmarked children to list of pairs, to be added to queue.
183  children.resize(0);
184  for (offset_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
185  gno_t child = edgeIds[ptr];
186  if (static_cast<Tpetra::global_size_t>(invPerm[child]) == INVALID){
187  // Not visited yet; add child to list of pairs.
188  std::pair<gno_t,offset_t> newchild;
189  newchild.first = child;
190  newchild.second = offsets[child+1] - offsets[child];
191  children.push_back(newchild);
192  }
193  }
194  // Sort children by increasing degree
195  // TODO: If edge weights, sort children by decreasing weight,
197  zort.sort(children);
198 
199  typename Teuchos::Array<std::pair<gno_t,offset_t> >::iterator it = children.begin ();
200  for ( ; it != children.end(); ++it){
201  // Push children on the queue in sorted order.
202  gno_t child = it->first;
203  invPerm[child] = count++; // Label as we push on Q
204  tmpPerm[invPerm[child]] = child;
205  Q.push(child);
206  //cout << "Debug: invPerm[" << child << "] = " << count << endl;
207  }
208  }
209  }
210 
211  // Old row tmpPerm[i] is now the new row i.
212 
213  // Reverse labels for RCM.
214  bool reverse = true; // TODO: Make parameter
215  if (reverse) {
216  lno_t temp;
217  for (size_t i=0; i < nVtx/2; ++i) {
218  // This effectively does the work of two loops:
219  // 1) for (i=1; i< nVtx/2; ++i)
220  // swap of tmpPerm[i] and tmpPerm[nVtx-1-i]
221  // 2) for (i=0; i < nVtx; ++i)
222  // invPerm[tmpPerm[i]] = i;
223  temp = tmpPerm[i];
224  invPerm[tmpPerm[nVtx-1-i]] = i;
225  invPerm[temp] = nVtx-1-i;
226  }
227 
228  }
229 
230  solution->setHaveInverse(true);
231  return ierr;
232  }
233 
234  private:
235  // Find a smallest degree vertex in component containing v
236  gno_t findSmallestDegree(
237  gno_t v,
238  lno_t nVtx,
239  ArrayView<const gno_t> edgeIds,
240  ArrayView<const offset_t> offsets)
241  {
242  std::queue<gno_t> Q;
243  Teuchos::Array<bool> mark(nVtx);
244 
245  // Do BFS and compute smallest degree as we go
246  offset_t smallestDegree = nVtx;
247  gno_t smallestVertex = 0;
248 
249  // Clear mark array - nothing marked yet
250  for (int i=0; i<nVtx; i++)
251  mark[i] = false;
252 
253  // Start from v
254  Q.push(v);
255  while (Q.size()){
256  // Get first vertex from the queue
257  v = Q.front();
258  Q.pop();
259  // Check degree of v
260  offset_t deg = offsets[v+1] - offsets[v];
261  if (deg < smallestDegree){
262  smallestDegree = deg;
263  smallestVertex = v;
264  }
265  // Add unmarked children to queue
266  for (offset_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
267  gno_t child = edgeIds[ptr];
268  if (!mark[child]){
269  mark[child] = true;
270  Q.push(child);
271  }
272  }
273  }
274  return smallestVertex;
275  }
276 
277  // Find a pseudoperipheral vertex in component containing v
278  gno_t findPseudoPeripheral(
279  gno_t v,
280  lno_t nVtx,
281  ArrayView<const gno_t> edgeIds,
282  ArrayView<const offset_t> offsets)
283  {
284  std::queue<gno_t> Q;
285  Teuchos::Array<bool> mark(nVtx);
286 
287  // Do BFS a couple times, pick vertex last visited (furthest away)
288  const int numBFS = 2;
289  for (int bfs=0; bfs<numBFS; bfs++){
290  // Clear mark array - nothing marked yet
291  for (int i=0; i<nVtx; i++)
292  mark[i] = false;
293  // Start from v
294  Q.push(v);
295  while (Q.size()){
296  // Get first vertex from the queue
297  v = Q.front();
298  Q.pop();
299  // Add unmarked children to queue
300  for (offset_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
301  gno_t child = edgeIds[ptr];
302  if (!mark[child]){
303  mark[child] = true;
304  Q.push(child);
305  }
306  }
307  }
308  }
309  return v;
310  }
311 
312 };
313 }
314 #endif
#define HELLO
#define INVALID(STR)
int localOrder(const RCP< LocalOrderingSolution< lno_t > > &solution)
Ordering method.
AlgRCM(const RCP< const typename Adapter::base_adapter_t > &adapter__, const RCP< Teuchos::ParameterList > &pl__, const RCP< const Teuchos::Comm< int > > &comm__, RCP< const Environment > &env__, const modelFlag_t &graphFlags__)
std::bitset< NUM_MODEL_FLAGS > modelFlag_t
map_t::global_ordinal_type gno_t
Definition: mapRemotes.cpp:18
Defines the OrderingSolution class.
Adapter::offset_t offset_t
tuple root
Definition: validXML.py:24
Algorithm defines the base class for all algorithms.
map_t::local_ordinal_type lno_t
Definition: mapRemotes.cpp:17
int globalOrder(const RCP< GlobalOrderingSolution< gno_t > > &)
Ordering method.
Sort vector of pairs (key, value) by value.
Adapter::scalar_t scalar_t
GraphModel defines the interface required for graph models.
Tpetra::global_size_t global_size_t
Defines the GraphModel interface.
Adapter::lno_t lno_t
void sort(Teuchos::Array< std::pair< key_t, value_t > > &listofPairs, bool inc=true)
Adapter::gno_t gno_t