EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_SymmRCM_CrsGraph.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
43#ifdef HAVE_EXPERIMENTAL
44#ifdef HAVE_GRAPH_REORDERINGS
45
47
49
50#include <set>
51
52#include <Epetra_Util.h>
53#include <Epetra_CrsGraph.h>
54#include <Epetra_Map.h>
55#include <Epetra_Import.h>
56
57namespace EpetraExt {
58
61{
62 if( newObj_ ) delete newObj_;
63
64 if( RCMColMap_ != RCMMap_ ) delete RCMColMap_;
65 if( RCMMap_ ) delete RCMMap_;
66}
67
71{
72 origObj_ = &orig;
73
74 int err;
75
76 //Generate Local Transpose Graph
77 CrsGraph_Transpose transposeTransform;
78 Epetra_CrsGraph & trans = transposeTransform( orig );
79
80 //Generate Local Symmetric Adj. List
81 //Find Min Degree Node While at it
82 int NumNodes = orig.NumMyRows();
83 int * LocalRow;
84 int * LocalRowTrans;
85 int RowSize, RowSizeTrans;
86 std::vector< std::vector<int> > AdjList( NumNodes );
87 int MinDegree = NumNodes;
88 int MinDegreeNode;
89 for( int i = 0; i < NumNodes; ++i )
90 {
91 orig.ExtractMyRowView( i, RowSize, LocalRow );
92 trans.ExtractMyRowView( i, RowSizeTrans, LocalRowTrans );
93
94 std::set<int> adjSet;
95 for( int j = 0; j < RowSize; ++j )
96 if( LocalRow[j] < NumNodes ) adjSet.insert( LocalRow[j] );
97 for( int j = 0; j < RowSizeTrans; ++j )
98 if( LocalRowTrans[j] < NumNodes ) adjSet.insert( LocalRowTrans[j] );
99
100 std::set<int>::iterator iterS = adjSet.begin();
101 std::set<int>::iterator endS = adjSet.end();
102 AdjList[i].resize( adjSet.size() );
103 for( int j = 0; iterS != endS; ++iterS, ++j ) AdjList[i][j] = *iterS;
104
105 if( AdjList[i].size() < MinDegree )
106 {
107 MinDegree = AdjList[i].size();
108 MinDegreeNode = i;
109 }
110 }
111
112 BFT * BestBFT;
113 bool TooWide;
114
115 //std::cout << "SymmRCM::bruteForce_ : " << bruteForce_ << std::endl;
116
117 if( bruteForce_ )
118 {
119 int bestWidth = NumNodes;
120 int bestDepth = 0;
121
122 for( int i = 0; i < NumNodes; ++i )
123 {
124 BFT * TestBFT = new BFT( AdjList, i, NumNodes, TooWide );
125 if( TestBFT->Depth() > bestDepth ||
126 ( TestBFT->Depth() == bestDepth && TestBFT->Width() < bestWidth ) )
127 {
128 BestBFT = TestBFT;
129 bestDepth = TestBFT->Depth();
130 bestWidth = TestBFT->Width();
131 }
132 else
133 delete TestBFT;
134 }
135 }
136 else
137 {
138 //Construct BFT for first
139 BestBFT = new BFT( AdjList, MinDegreeNode, NumNodes, TooWide );
140
141 int MinWidth = BestBFT->Width();
142 int BestWidth = MinWidth;
143 int Diameter = BestBFT->Depth();
144 std::vector<int> Leaves;
145 BestBFT->NonNeighborLeaves( Leaves, AdjList, testLeafWidth_ );
146
147 bool DeeperFound;
148 bool NarrowerFound;
149
150 bool Finished = false;
151
152 while( !Finished )
153 {
154 DeeperFound = false;
155 NarrowerFound = false;
156
157 for( int i = 0; i < Leaves.size(); ++i )
158 {
159
160 BFT * TestBFT = new BFT( AdjList, Leaves[i], MinWidth, TooWide );
161
162 if( TooWide )
163 delete TestBFT;
164 else
165 {
166 if( TestBFT->Width() < MinWidth ) MinWidth = TestBFT->Width();
167
168 if( TestBFT->Depth() > Diameter )
169 {
170 delete BestBFT;
171 Diameter = TestBFT->Depth();
172 BestWidth = TestBFT->Width();
173 BestBFT = TestBFT;
174 DeeperFound = true;
175 NarrowerFound = false;
176 }
177 else if( (TestBFT->Depth()==Diameter) && (TestBFT->Width()<BestWidth) )
178 {
179 delete BestBFT;
180 BestWidth = TestBFT->Width();
181 BestBFT = TestBFT;
182 NarrowerFound = true;
183 }
184 else delete TestBFT;
185 }
186 }
187
188 if( DeeperFound )
189 BestBFT->NonNeighborLeaves( Leaves, AdjList, testLeafWidth_ );
190 else if( NarrowerFound )
191 Finished = true;
192 else Finished = true;
193 }
194 }
195
196 //std::cout << "\nSymmRCM:\n";
197 //std::cout << "----------------------------\n";
198 //std::cout << " Depth: " << BestBFT->Depth() << std::endl;
199 //std::cout << " Width: " << BestBFT->Width() << std::endl;
200 //std::cout << "----------------------------\n\n";
201
202 std::vector<int> RCM;
203 BestBFT->ReverseVector( RCM );
204 for( int i = 0; i < NumNodes; ++i )
205 RCM[i] = orig.RowMap().GID( RCM[i] );
206
207 //Generate New Row Map
208 RCMMap_ = new Epetra_Map( orig.RowMap().NumGlobalElements(),
209 NumNodes,
210 &RCM[0],
211 orig.RowMap().IndexBase(),
212 orig.RowMap().Comm() );
213
214 //Generate New Col Map
215 if( RCMMap_->DistributedGlobal() )
216 {
217 std::vector<int> colIndices = RCM;
218 const Epetra_BlockMap & origColMap = orig.ColMap();
219
220 if( origColMap.NumMyElements() > RCMMap_->NumMyElements() )
221 {
222 for( int i = RCMMap_->NumMyElements(); i < origColMap.NumMyElements(); ++i )
223 colIndices.push_back( origColMap.GID(i) );
224 }
225
226 RCMColMap_ = new Epetra_Map( orig.ColMap().NumGlobalElements(),
227 colIndices.size(),
228 &colIndices[0],
229 orig.ColMap().IndexBase(),
230 orig.ColMap().Comm() );
231 }
232 else
233 RCMColMap_ = RCMMap_;
234
235 //Create New Graph
236 Epetra_Import Importer( *RCMMap_, orig.RowMap() );
237 Epetra_CrsGraph * RCMGraph = new Epetra_CrsGraph( Copy, *RCMMap_, *RCMColMap_, 0 );
238 RCMGraph->Import( orig, Importer, Insert );
239 RCMGraph->FillComplete();
240
241/*
242 std::cout << "origGraph\n";
243 std::cout << orig;
244 std::cout << "RCMGraph\n";
245 std::cout << *RCMGraph;
246*/
247
248 newObj_ = RCMGraph;
249
250 return *RCMGraph;
251}
252
253CrsGraph_SymmRCM::BFT::
254BFT( const std::vector< std::vector<int> > & adjlist,
255 int root,
256 int max_width,
257 bool & failed )
258: width_(0),
259 depth_(0),
260 nodes_(adjlist.size()),
261 failed_(false)
262{
263 std::set<int> touchedNodes;
264
265 //setup level 0 of traversal
266 levelSets_.push_back( std::vector<int>(1) );
267 levelSets_[0][0] = root;
268 ++depth_;
269
270 //start set of touched nodes
271 touchedNodes.insert( root );
272
273 while( touchedNodes.size() < nodes_ )
274 {
275 //start new level set
276 levelSets_.push_back( std::vector<int>() );
277 ++depth_;
278
279 for( int i = 0; i < levelSets_[depth_-2].size(); ++i )
280 {
281 int currNode = levelSets_[depth_-2][i];
282 int adjSize = adjlist[currNode].size();
283 for( int j = 0; j < adjSize; ++j )
284 {
285 // add nodes to current level set when new
286 int currAdj = adjlist[currNode][j];
287 if( !touchedNodes.count( currAdj ) )
288 {
289 levelSets_[depth_-1].push_back( currAdj );
290 touchedNodes.insert( currAdj );
291 }
292 }
293 }
294
295 int currWidth = levelSets_[depth_-1].size();
296
297 if( currWidth ) //sort adj nodes by degree
298 {
299 if( currWidth > width_ ) width_ = currWidth;
300
301 //fail if width is greater than allowed
302 if( width_ > max_width )
303 {
304 failed_ = true;
305 failed = failed_;
306 return;
307 }
308
309 //Increasing Order By Degree
310 std::vector<int> degrees( currWidth );
311 for( int i = 0; i < currWidth; ++i )
312 degrees[i] = adjlist[ levelSets_[depth_-1][i] ].size();
313 int ** indices = new int*[1];
314 indices[0] = &(levelSets_[depth_-1][0]);
315 Epetra_Util().Sort( true, currWidth, &degrees[0], 0, 0, 1, indices );
316 }
317 else //it is a disconnected graph
318 {
319 //start again from minimum degree node of those remaining
320 bool found = false;
321 int minDegree = nodes_;
322 int minDegreeNode;
323 for( int i = 0; i < nodes_; ++i )
324 {
325 if( !touchedNodes.count( i ) && adjlist[i].size() < minDegree )
326 {
327 minDegree = adjlist[i].size();
328 minDegreeNode = i;
329 found = true;
330 }
331 }
332
333 if( found )
334 {
335 touchedNodes.insert( minDegreeNode );
336 levelSets_[depth_-1].push_back( minDegreeNode );
337 }
338 else
339 {
340 --depth_;
341 failed_ = true;
342 failed = failed_;
343 return;
344 }
345
346 }
347
348 }
349
350/*
351 std::cout << "BFT<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
352 std::cout << "Width: " << width_ << std::endl;
353 std::cout << "Depth: " << depth_ << std::endl;
354 std::cout << "Adj List: " << nodes_ << std::endl;
355 for( int i = 0; i < nodes_; ++i )
356 {
357 std::cout << i << "\t";
358 for( int j = 0; j < adjlist[i].size(); ++j )
359 std::cout << adjlist[i][j] << " ";
360 std::cout << std::endl;
361 }
362 std::cout << "Level Sets: " << depth_ << std::endl;
363 for( int i = 0; i < depth_; ++i )
364 {
365 std::cout << i << "\t";
366 for( int j = 0; j < levelSets_[i].size(); ++j )
367 std::cout << levelSets_[i][j] << " ";
368 std::cout << std::endl;
369 }
370 std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
371*/
372
373 failed = failed_;
374}
375
376void
377CrsGraph_SymmRCM::BFT::
378NonNeighborLeaves( std::vector<int> & leaves,
379 const std::vector< std::vector<int> > & adjlist,
380 int count )
381{
382 assert( (depth_>0) && (failed_==false) );
383
384 leaves.clear();
385 int leafWidth = levelSets_[depth_-1].size();
386 std::set<int> adjSeen;
387 for( int i = 0; i < leafWidth; ++i )
388 {
389 int currLeaf = levelSets_[depth_-1][i];
390 if( !adjSeen.count( currLeaf ) )
391 {
392 leaves.push_back( currLeaf );
393 for( int j = 0; j < adjlist[currLeaf].size(); ++j )
394 adjSeen.insert( adjlist[currLeaf][j] );
395 }
396 if( leaves.size() == count ) i = leafWidth;
397 }
398}
399
400void
401CrsGraph_SymmRCM::BFT::
402ReverseVector( std::vector<int> & ordered )
403{
404 assert( (depth_>0) && (failed_==false) );
405
406 ordered.resize( nodes_ );
407 int loc = 0;
408 for( int i = 0; i < depth_; ++i )
409 {
410 int currLevel = depth_ - (i+1);
411 int currWidth = levelSets_[currLevel].size();
412 for( int j = 0; j < currWidth; ++j )
413 ordered[loc++] = levelSets_[currLevel][currWidth-j-1];
414 }
415}
416
417} //namespace EpetraExt
418#endif //HAVE_GRAPH_REORDERINGS
419#endif //HAVE_EXPERIMENTAL
bool DistributedGlobal() const
int NumMyElements() const
static void EPETRA_LIB_DLL_EXPORT Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
NewTypeRef operator()(OriginalTypeRef orig)
Transformation Operator.