EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_SymmRCM_CrsGraph.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include "EpetraExt_ConfigDefs.h"
00043 #ifdef HAVE_EXPERIMENTAL
00044 #ifdef HAVE_GRAPH_REORDERINGS
00045 
00046 #include <EpetraExt_SymmRCM_CrsGraph.h>
00047 
00048 #include <EpetraExt_Transpose_CrsGraph.h>
00049 
00050 #include <set>
00051 
00052 #include <Epetra_Util.h>
00053 #include <Epetra_CrsGraph.h>
00054 #include <Epetra_Map.h>
00055 #include <Epetra_Import.h>
00056 
00057 namespace EpetraExt {
00058 
00059 CrsGraph_SymmRCM::
00060 ~CrsGraph_SymmRCM()
00061 {
00062   if( newObj_ ) delete newObj_;
00063 
00064   if( RCMColMap_ != RCMMap_ ) delete RCMColMap_;
00065   if( RCMMap_ ) delete RCMMap_;
00066 }
00067 
00068 CrsGraph_SymmRCM::NewTypeRef
00069 CrsGraph_SymmRCM::
00070 operator()( CrsGraph_SymmRCM::OriginalTypeRef orig )
00071 {
00072   origObj_ = &orig;
00073 
00074   int err;
00075 
00076   //Generate Local Transpose Graph
00077   CrsGraph_Transpose transposeTransform;
00078   Epetra_CrsGraph & trans = transposeTransform( orig );
00079 
00080   //Generate Local Symmetric Adj. List
00081   //Find Min Degree Node While at it
00082   int NumNodes = orig.NumMyRows();
00083   int * LocalRow;
00084   int * LocalRowTrans;
00085   int RowSize, RowSizeTrans;
00086   std::vector< std::vector<int> > AdjList( NumNodes );
00087   int MinDegree = NumNodes;
00088   int MinDegreeNode;
00089   for( int i = 0; i < NumNodes; ++i )
00090   {
00091     orig.ExtractMyRowView( i, RowSize, LocalRow );
00092     trans.ExtractMyRowView( i, RowSizeTrans, LocalRowTrans );
00093 
00094     std::set<int> adjSet;
00095     for( int j = 0; j < RowSize; ++j )
00096      if( LocalRow[j] < NumNodes ) adjSet.insert( LocalRow[j] );
00097     for( int j = 0; j < RowSizeTrans; ++j )
00098      if( LocalRowTrans[j] < NumNodes ) adjSet.insert( LocalRowTrans[j] );
00099 
00100     std::set<int>::iterator iterS = adjSet.begin();
00101     std::set<int>::iterator endS = adjSet.end();
00102     AdjList[i].resize( adjSet.size() );
00103     for( int j = 0; iterS != endS; ++iterS, ++j ) AdjList[i][j] = *iterS;
00104     
00105     if( AdjList[i].size() < MinDegree )
00106     {
00107       MinDegree = AdjList[i].size();
00108       MinDegreeNode = i;
00109     }
00110   }
00111 
00112   BFT * BestBFT;
00113   bool TooWide;
00114 
00115   //std::cout << "SymmRCM::bruteForce_ : " << bruteForce_ << std::endl;
00116 
00117   if( bruteForce_ )
00118   {
00119     int bestWidth = NumNodes;
00120     int bestDepth = 0;
00121     
00122     for( int i = 0; i < NumNodes; ++i )
00123     {
00124       BFT * TestBFT = new BFT( AdjList, i, NumNodes, TooWide );
00125       if( TestBFT->Depth() > bestDepth ||
00126           ( TestBFT->Depth() == bestDepth && TestBFT->Width() < bestWidth ) )
00127       {
00128         BestBFT = TestBFT;
00129         bestDepth = TestBFT->Depth();
00130         bestWidth = TestBFT->Width();
00131       }
00132       else
00133         delete TestBFT;
00134     }
00135   }
00136   else
00137   {
00138     //Construct BFT for first
00139     BestBFT = new BFT( AdjList, MinDegreeNode, NumNodes, TooWide );
00140 
00141     int MinWidth = BestBFT->Width();
00142     int BestWidth = MinWidth;
00143     int Diameter = BestBFT->Depth();
00144     std::vector<int> Leaves;
00145     BestBFT->NonNeighborLeaves( Leaves, AdjList, testLeafWidth_ );
00146 
00147     bool DeeperFound;
00148     bool NarrowerFound;
00149   
00150     bool Finished = false;
00151 
00152     while( !Finished )
00153     {
00154       DeeperFound = false;
00155       NarrowerFound = false;
00156 
00157       for( int i = 0; i < Leaves.size(); ++i )
00158       {
00159 
00160         BFT * TestBFT = new BFT( AdjList, Leaves[i], MinWidth, TooWide );
00161 
00162         if( TooWide )
00163           delete TestBFT;
00164         else
00165         {
00166           if( TestBFT->Width() < MinWidth ) MinWidth = TestBFT->Width();
00167 
00168           if( TestBFT->Depth() > Diameter )
00169           {
00170             delete BestBFT;
00171             Diameter = TestBFT->Depth();
00172             BestWidth = TestBFT->Width();
00173             BestBFT = TestBFT;
00174             DeeperFound = true;
00175             NarrowerFound = false;
00176           }
00177           else if( (TestBFT->Depth()==Diameter) && (TestBFT->Width()<BestWidth) )
00178           {
00179             delete BestBFT;
00180             BestWidth = TestBFT->Width();
00181             BestBFT = TestBFT;
00182             NarrowerFound = true;
00183           }
00184           else delete TestBFT;
00185         }
00186       }
00187 
00188       if( DeeperFound )
00189         BestBFT->NonNeighborLeaves( Leaves, AdjList, testLeafWidth_ );
00190       else if( NarrowerFound )
00191         Finished = true;
00192       else Finished = true;
00193     }
00194   }
00195 
00196   //std::cout << "\nSymmRCM:\n";
00197   //std::cout << "----------------------------\n";
00198   //std::cout << " Depth: " << BestBFT->Depth() << std::endl;
00199   //std::cout << " Width: " << BestBFT->Width() << std::endl;
00200   //std::cout << "----------------------------\n\n";
00201 
00202   std::vector<int> RCM;
00203   BestBFT->ReverseVector( RCM );
00204   for( int i = 0; i < NumNodes; ++i )
00205     RCM[i] = orig.RowMap().GID( RCM[i] );
00206 
00207   //Generate New Row Map
00208   RCMMap_ = new Epetra_Map( orig.RowMap().NumGlobalElements(),
00209                                         NumNodes,
00210                                         &RCM[0],
00211                                         orig.RowMap().IndexBase(),
00212                                         orig.RowMap().Comm() );
00213 
00214   //Generate New Col Map
00215   if( RCMMap_->DistributedGlobal() )
00216   {
00217     std::vector<int> colIndices = RCM;
00218     const Epetra_BlockMap & origColMap = orig.ColMap();
00219 
00220     if( origColMap.NumMyElements() > RCMMap_->NumMyElements() )
00221     {
00222       for( int i = RCMMap_->NumMyElements(); i < origColMap.NumMyElements(); ++i )
00223         colIndices.push_back( origColMap.GID(i) );
00224     }
00225 
00226     RCMColMap_ = new Epetra_Map( orig.ColMap().NumGlobalElements(),
00227                                  colIndices.size(),
00228                                  &colIndices[0],
00229                                  orig.ColMap().IndexBase(),
00230                                  orig.ColMap().Comm() );
00231   }
00232   else
00233     RCMColMap_ = RCMMap_;
00234 
00235   //Create New Graph
00236   Epetra_Import Importer( *RCMMap_, orig.RowMap() );
00237   Epetra_CrsGraph * RCMGraph = new Epetra_CrsGraph( Copy, *RCMMap_, *RCMColMap_, 0 );
00238   RCMGraph->Import( orig, Importer, Insert );
00239   RCMGraph->FillComplete();
00240 
00241 /*
00242   std::cout << "origGraph\n";
00243   std::cout << orig;
00244   std::cout << "RCMGraph\n";
00245   std::cout << *RCMGraph;
00246 */
00247 
00248   newObj_ = RCMGraph;
00249   
00250   return *RCMGraph;
00251 }
00252 
00253 CrsGraph_SymmRCM::BFT::
00254 BFT( const std::vector< std::vector<int> > & adjlist,
00255      int root,
00256      int max_width,
00257      bool & failed )
00258 : width_(0),
00259   depth_(0),
00260   nodes_(adjlist.size()),
00261   failed_(false)
00262 {
00263   std::set<int> touchedNodes;
00264 
00265   //setup level 0 of traversal
00266   levelSets_.push_back( std::vector<int>(1) );
00267   levelSets_[0][0] = root;
00268   ++depth_;
00269 
00270   //start set of touched nodes
00271   touchedNodes.insert( root );
00272 
00273   while( touchedNodes.size() < nodes_ )
00274   {
00275     //start new level set
00276     levelSets_.push_back( std::vector<int>() );
00277     ++depth_;
00278 
00279     for( int i = 0; i < levelSets_[depth_-2].size(); ++i )
00280     {
00281       int currNode = levelSets_[depth_-2][i];
00282       int adjSize  = adjlist[currNode].size();
00283       for( int j = 0; j < adjSize; ++j )
00284       {
00285         // add nodes to current level set when new
00286         int currAdj = adjlist[currNode][j];
00287         if( !touchedNodes.count( currAdj ) )
00288         {
00289           levelSets_[depth_-1].push_back( currAdj );
00290           touchedNodes.insert( currAdj );
00291         }
00292       }
00293     }
00294 
00295     int currWidth = levelSets_[depth_-1].size();
00296 
00297     if( currWidth ) //sort adj nodes by degree
00298     {
00299       if( currWidth > width_ ) width_ = currWidth;
00300 
00301       //fail if width is greater than allowed
00302       if( width_ > max_width )
00303       {
00304         failed_ = true;
00305         failed = failed_;
00306         return;
00307       }
00308 
00309       //Increasing Order By Degree
00310       std::vector<int> degrees( currWidth );
00311       for( int i = 0; i < currWidth; ++i )
00312         degrees[i] = adjlist[ levelSets_[depth_-1][i] ].size();
00313       int ** indices = new int*[1];
00314       indices[0] = &(levelSets_[depth_-1][0]);
00315       Epetra_Util().Sort( true, currWidth, &degrees[0], 0, 0, 1, indices );
00316     }
00317     else //it is a disconnected graph
00318     {
00319       //start again from minimum degree node of those remaining
00320       bool found = false;
00321       int minDegree = nodes_;
00322       int minDegreeNode;
00323       for( int i = 0; i < nodes_; ++i )
00324       {
00325         if( !touchedNodes.count( i ) && adjlist[i].size() < minDegree )
00326         {
00327           minDegree = adjlist[i].size();
00328           minDegreeNode = i;
00329           found = true;
00330         }
00331       }
00332 
00333       if( found )
00334       {
00335         touchedNodes.insert( minDegreeNode );
00336         levelSets_[depth_-1].push_back( minDegreeNode );
00337       }
00338       else
00339       {
00340         --depth_;
00341         failed_ = true;
00342         failed = failed_;
00343         return;
00344       }
00345 
00346     }
00347 
00348   }
00349 
00350 /*
00351   std::cout << "BFT<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
00352   std::cout << "Width: " << width_ << std::endl;
00353   std::cout << "Depth: " << depth_ << std::endl;
00354   std::cout << "Adj List: " << nodes_ << std::endl;
00355   for( int i = 0; i < nodes_; ++i )
00356   {
00357     std::cout << i << "\t";
00358     for( int j = 0; j < adjlist[i].size(); ++j )
00359       std::cout << adjlist[i][j] << " ";
00360     std::cout << std::endl;
00361   }
00362   std::cout << "Level Sets: " << depth_ << std::endl;
00363   for( int i = 0; i < depth_; ++i )
00364   {
00365     std::cout << i << "\t";
00366     for( int j = 0; j < levelSets_[i].size(); ++j )
00367       std::cout << levelSets_[i][j] << " ";
00368     std::cout << std::endl;
00369   }
00370   std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
00371 */
00372 
00373   failed = failed_;
00374 }
00375 
00376 void
00377 CrsGraph_SymmRCM::BFT::
00378 NonNeighborLeaves( std::vector<int> & leaves,
00379                    const std::vector< std::vector<int> > & adjlist,
00380                    int count )
00381 {
00382   assert( (depth_>0) && (failed_==false) );
00383 
00384   leaves.clear();
00385   int leafWidth = levelSets_[depth_-1].size();
00386   std::set<int> adjSeen;
00387   for( int i = 0; i < leafWidth; ++i )
00388   {
00389     int currLeaf = levelSets_[depth_-1][i];
00390     if( !adjSeen.count( currLeaf ) )
00391     {
00392       leaves.push_back( currLeaf );
00393       for( int j = 0; j < adjlist[currLeaf].size(); ++j )
00394         adjSeen.insert( adjlist[currLeaf][j] );
00395     }
00396     if( leaves.size() == count ) i = leafWidth;
00397   }
00398 }
00399 
00400 void
00401 CrsGraph_SymmRCM::BFT::
00402 ReverseVector( std::vector<int> & ordered )
00403 {
00404   assert( (depth_>0) && (failed_==false) );
00405 
00406   ordered.resize( nodes_ );
00407   int loc = 0;
00408   for( int i = 0; i < depth_; ++i )
00409   {
00410     int currLevel = depth_ - (i+1);
00411     int currWidth = levelSets_[currLevel].size();
00412     for( int j = 0; j < currWidth; ++j )
00413       ordered[loc++] = levelSets_[currLevel][currWidth-j-1];
00414   }
00415 }
00416 
00417 } //namespace EpetraExt
00418 #endif //HAVE_GRAPH_REORDERINGS
00419 #endif //HAVE_EXPERIMENTAL
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines