EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
42 #include "EpetraExt_ConfigDefs.h"
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 
57 namespace 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
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 
254 BFT( 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 
376 void
378 NonNeighborLeaves( 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 
400 void
402 ReverseVector( 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
NewTypeRef operator()(OriginalTypeRef orig)
Transformation Operator.
std::vector< std::vector< int > > levelSets_
void NonNeighborLeaves(std::vector< int > &leaves, const std::vector< std::vector< int > > &adjlist, int count)
int NumMyElements() const
~CrsGraph_SymmRCM()
Destructor.
void ReverseVector(std::vector< int > &ordered)
int GID(int LID) const
BFT(const std::vector< std::vector< int > > &adjlist, int root, int max_width, bool &failed)
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
static void Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)