Zoltan2
Loading...
Searching...
No Matches
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>
51#include <Zoltan2_Sort.hpp>
52#include <queue>
53
54
58
59
60namespace Zoltan2{
61
62template <typename Adapter>
63class 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,
196 SortPairs<gno_t,offset_t> zort;
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 INVALID(STR)
Defines the GraphModel interface.
Defines the OrderingSolution class.
Sort vector of pairs (key, value) by value.
#define HELLO
Adapter::gno_t gno_t
Adapter::lno_t lno_t
Adapter::offset_t offset_t
int localOrder(const RCP< LocalOrderingSolution< lno_t > > &solution)
Ordering method.
int globalOrder(const RCP< GlobalOrderingSolution< gno_t > > &)
Ordering method.
Adapter::scalar_t scalar_t
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__)
Algorithm defines the base class for all algorithms.
map_t::local_ordinal_type lno_t
map_t::global_ordinal_type gno_t
Created by mbenlioglu on Aug 31, 2020.
std::bitset< NUM_MODEL_FLAGS > modelFlag_t