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 (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "EpetraExt_ConfigDefs.h"
00030 #ifdef HAVE_EXPERIMENTAL
00031 #ifdef HAVE_GRAPH_REORDERINGS
00032 
00033 #include <EpetraExt_SymmRCM_CrsGraph.h>
00034 
00035 #include <EpetraExt_Transpose_CrsGraph.h>
00036 
00037 #include <set>
00038 
00039 #include <Epetra_Util.h>
00040 #include <Epetra_CrsGraph.h>
00041 #include <Epetra_Map.h>
00042 #include <Epetra_Import.h>
00043 
00044 namespace EpetraExt {
00045 
00046 CrsGraph_SymmRCM::
00047 ~CrsGraph_SymmRCM()
00048 {
00049   if( newObj_ ) delete newObj_;
00050 
00051   if( RCMColMap_ != RCMMap_ ) delete RCMColMap_;
00052   if( RCMMap_ ) delete RCMMap_;
00053 }
00054 
00055 CrsGraph_SymmRCM::NewTypeRef
00056 CrsGraph_SymmRCM::
00057 operator()( CrsGraph_SymmRCM::OriginalTypeRef orig )
00058 {
00059   origObj_ = &orig;
00060 
00061   int err;
00062 
00063   //Generate Local Transpose Graph
00064   CrsGraph_Transpose transposeTransform;
00065   Epetra_CrsGraph & trans = transposeTransform( orig );
00066 
00067   //Generate Local Symmetric Adj. List
00068   //Find Min Degree Node While at it
00069   int NumNodes = orig.NumMyRows();
00070   int * LocalRow;
00071   int * LocalRowTrans;
00072   int RowSize, RowSizeTrans;
00073   std::vector< std::vector<int> > AdjList( NumNodes );
00074   int MinDegree = NumNodes;
00075   int MinDegreeNode;
00076   for( int i = 0; i < NumNodes; ++i )
00077   {
00078     orig.ExtractMyRowView( i, RowSize, LocalRow );
00079     trans.ExtractMyRowView( i, RowSizeTrans, LocalRowTrans );
00080 
00081     std::set<int> adjSet;
00082     for( int j = 0; j < RowSize; ++j )
00083      if( LocalRow[j] < NumNodes ) adjSet.insert( LocalRow[j] );
00084     for( int j = 0; j < RowSizeTrans; ++j )
00085      if( LocalRowTrans[j] < NumNodes ) adjSet.insert( LocalRowTrans[j] );
00086 
00087     std::set<int>::iterator iterS = adjSet.begin();
00088     std::set<int>::iterator endS = adjSet.end();
00089     AdjList[i].resize( adjSet.size() );
00090     for( int j = 0; iterS != endS; ++iterS, ++j ) AdjList[i][j] = *iterS;
00091     
00092     if( AdjList[i].size() < MinDegree )
00093     {
00094       MinDegree = AdjList[i].size();
00095       MinDegreeNode = i;
00096     }
00097   }
00098 
00099   BFT * BestBFT;
00100   bool TooWide;
00101 
00102   //std::cout << "SymmRCM::bruteForce_ : " << bruteForce_ << std::endl;
00103 
00104   if( bruteForce_ )
00105   {
00106     int bestWidth = NumNodes;
00107     int bestDepth = 0;
00108     
00109     for( int i = 0; i < NumNodes; ++i )
00110     {
00111       BFT * TestBFT = new BFT( AdjList, i, NumNodes, TooWide );
00112       if( TestBFT->Depth() > bestDepth ||
00113           ( TestBFT->Depth() == bestDepth && TestBFT->Width() < bestWidth ) )
00114       {
00115         BestBFT = TestBFT;
00116         bestDepth = TestBFT->Depth();
00117         bestWidth = TestBFT->Width();
00118       }
00119       else
00120         delete TestBFT;
00121     }
00122   }
00123   else
00124   {
00125     //Construct BFT for first
00126     BestBFT = new BFT( AdjList, MinDegreeNode, NumNodes, TooWide );
00127 
00128     int MinWidth = BestBFT->Width();
00129     int BestWidth = MinWidth;
00130     int Diameter = BestBFT->Depth();
00131     std::vector<int> Leaves;
00132     BestBFT->NonNeighborLeaves( Leaves, AdjList, testLeafWidth_ );
00133 
00134     bool DeeperFound;
00135     bool NarrowerFound;
00136   
00137     bool Finished = false;
00138 
00139     while( !Finished )
00140     {
00141       DeeperFound = false;
00142       NarrowerFound = false;
00143 
00144       for( int i = 0; i < Leaves.size(); ++i )
00145       {
00146 
00147         BFT * TestBFT = new BFT( AdjList, Leaves[i], MinWidth, TooWide );
00148 
00149         if( TooWide )
00150           delete TestBFT;
00151         else
00152         {
00153           if( TestBFT->Width() < MinWidth ) MinWidth = TestBFT->Width();
00154 
00155           if( TestBFT->Depth() > Diameter )
00156           {
00157             delete BestBFT;
00158             Diameter = TestBFT->Depth();
00159             BestWidth = TestBFT->Width();
00160             BestBFT = TestBFT;
00161             DeeperFound = true;
00162             NarrowerFound = false;
00163           }
00164           else if( (TestBFT->Depth()==Diameter) && (TestBFT->Width()<BestWidth) )
00165           {
00166             delete BestBFT;
00167             BestWidth = TestBFT->Width();
00168             BestBFT = TestBFT;
00169             NarrowerFound = true;
00170           }
00171           else delete TestBFT;
00172         }
00173       }
00174 
00175       if( DeeperFound )
00176         BestBFT->NonNeighborLeaves( Leaves, AdjList, testLeafWidth_ );
00177       else if( NarrowerFound )
00178         Finished = true;
00179       else Finished = true;
00180     }
00181   }
00182 
00183   //std::cout << "\nSymmRCM:\n";
00184   //std::cout << "----------------------------\n";
00185   //std::cout << " Depth: " << BestBFT->Depth() << std::endl;
00186   //std::cout << " Width: " << BestBFT->Width() << std::endl;
00187   //std::cout << "----------------------------\n\n";
00188 
00189   std::vector<int> RCM;
00190   BestBFT->ReverseVector( RCM );
00191   for( int i = 0; i < NumNodes; ++i )
00192     RCM[i] = orig.RowMap().GID( RCM[i] );
00193 
00194   //Generate New Row Map
00195   RCMMap_ = new Epetra_Map( orig.RowMap().NumGlobalElements(),
00196                                         NumNodes,
00197                                         &RCM[0],
00198                                         orig.RowMap().IndexBase(),
00199                                         orig.RowMap().Comm() );
00200 
00201   //Generate New Col Map
00202   if( RCMMap_->DistributedGlobal() )
00203   {
00204     std::vector<int> colIndices = RCM;
00205     const Epetra_BlockMap & origColMap = orig.ColMap();
00206 
00207     if( origColMap.NumMyElements() > RCMMap_->NumMyElements() )
00208     {
00209       for( int i = RCMMap_->NumMyElements(); i < origColMap.NumMyElements(); ++i )
00210         colIndices.push_back( origColMap.GID(i) );
00211     }
00212 
00213     RCMColMap_ = new Epetra_Map( orig.ColMap().NumGlobalElements(),
00214                                  colIndices.size(),
00215                                  &colIndices[0],
00216                                  orig.ColMap().IndexBase(),
00217                                  orig.ColMap().Comm() );
00218   }
00219   else
00220     RCMColMap_ = RCMMap_;
00221 
00222   //Create New Graph
00223   Epetra_Import Importer( *RCMMap_, orig.RowMap() );
00224   Epetra_CrsGraph * RCMGraph = new Epetra_CrsGraph( Copy, *RCMMap_, *RCMColMap_, 0 );
00225   RCMGraph->Import( orig, Importer, Insert );
00226   RCMGraph->FillComplete();
00227 
00228 /*
00229   std::cout << "origGraph\n";
00230   std::cout << orig;
00231   std::cout << "RCMGraph\n";
00232   std::cout << *RCMGraph;
00233 */
00234 
00235   newObj_ = RCMGraph;
00236   
00237   return *RCMGraph;
00238 }
00239 
00240 CrsGraph_SymmRCM::BFT::
00241 BFT( const std::vector< std::vector<int> > & adjlist,
00242      int root,
00243      int max_width,
00244      bool & failed )
00245 : width_(0),
00246   depth_(0),
00247   nodes_(adjlist.size()),
00248   failed_(false)
00249 {
00250   std::set<int> touchedNodes;
00251 
00252   //setup level 0 of traversal
00253   levelSets_.push_back( std::vector<int>(1) );
00254   levelSets_[0][0] = root;
00255   ++depth_;
00256 
00257   //start set of touched nodes
00258   touchedNodes.insert( root );
00259 
00260   while( touchedNodes.size() < nodes_ )
00261   {
00262     //start new level set
00263     levelSets_.push_back( std::vector<int>() );
00264     ++depth_;
00265 
00266     for( int i = 0; i < levelSets_[depth_-2].size(); ++i )
00267     {
00268       int currNode = levelSets_[depth_-2][i];
00269       int adjSize  = adjlist[currNode].size();
00270       for( int j = 0; j < adjSize; ++j )
00271       {
00272         // add nodes to current level set when new
00273         int currAdj = adjlist[currNode][j];
00274         if( !touchedNodes.count( currAdj ) )
00275         {
00276           levelSets_[depth_-1].push_back( currAdj );
00277           touchedNodes.insert( currAdj );
00278         }
00279       }
00280     }
00281 
00282     int currWidth = levelSets_[depth_-1].size();
00283 
00284     if( currWidth ) //sort adj nodes by degree
00285     {
00286       if( currWidth > width_ ) width_ = currWidth;
00287 
00288       //fail if width is greater than allowed
00289       if( width_ > max_width )
00290       {
00291         failed_ = true;
00292         failed = failed_;
00293         return;
00294       }
00295 
00296       //Increasing Order By Degree
00297       std::vector<int> degrees( currWidth );
00298       for( int i = 0; i < currWidth; ++i )
00299         degrees[i] = adjlist[ levelSets_[depth_-1][i] ].size();
00300       int ** indices = new int*[1];
00301       indices[0] = &(levelSets_[depth_-1][0]);
00302       Epetra_Util().Sort( true, currWidth, &degrees[0], 0, 0, 1, indices );
00303     }
00304     else //it is a disconnected graph
00305     {
00306       //start again from minimum degree node of those remaining
00307       bool found = false;
00308       int minDegree = nodes_;
00309       int minDegreeNode;
00310       for( int i = 0; i < nodes_; ++i )
00311       {
00312         if( !touchedNodes.count( i ) && adjlist[i].size() < minDegree )
00313         {
00314           minDegree = adjlist[i].size();
00315           minDegreeNode = i;
00316           found = true;
00317         }
00318       }
00319 
00320       if( found )
00321       {
00322         touchedNodes.insert( minDegreeNode );
00323         levelSets_[depth_-1].push_back( minDegreeNode );
00324       }
00325       else
00326       {
00327         --depth_;
00328         failed_ = true;
00329         failed = failed_;
00330         return;
00331       }
00332 
00333     }
00334 
00335   }
00336 
00337 /*
00338   std::cout << "BFT<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
00339   std::cout << "Width: " << width_ << std::endl;
00340   std::cout << "Depth: " << depth_ << std::endl;
00341   std::cout << "Adj List: " << nodes_ << std::endl;
00342   for( int i = 0; i < nodes_; ++i )
00343   {
00344     std::cout << i << "\t";
00345     for( int j = 0; j < adjlist[i].size(); ++j )
00346       std::cout << adjlist[i][j] << " ";
00347     std::cout << std::endl;
00348   }
00349   std::cout << "Level Sets: " << depth_ << std::endl;
00350   for( int i = 0; i < depth_; ++i )
00351   {
00352     std::cout << i << "\t";
00353     for( int j = 0; j < levelSets_[i].size(); ++j )
00354       std::cout << levelSets_[i][j] << " ";
00355     std::cout << std::endl;
00356   }
00357   std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
00358 */
00359 
00360   failed = failed_;
00361 }
00362 
00363 void
00364 CrsGraph_SymmRCM::BFT::
00365 NonNeighborLeaves( std::vector<int> & leaves,
00366                    const std::vector< std::vector<int> > & adjlist,
00367                    int count )
00368 {
00369   assert( (depth_>0) && (failed_==false) );
00370 
00371   leaves.clear();
00372   int leafWidth = levelSets_[depth_-1].size();
00373   std::set<int> adjSeen;
00374   for( int i = 0; i < leafWidth; ++i )
00375   {
00376     int currLeaf = levelSets_[depth_-1][i];
00377     if( !adjSeen.count( currLeaf ) )
00378     {
00379       leaves.push_back( currLeaf );
00380       for( int j = 0; j < adjlist[currLeaf].size(); ++j )
00381         adjSeen.insert( adjlist[currLeaf][j] );
00382     }
00383     if( leaves.size() == count ) i = leafWidth;
00384   }
00385 }
00386 
00387 void
00388 CrsGraph_SymmRCM::BFT::
00389 ReverseVector( std::vector<int> & ordered )
00390 {
00391   assert( (depth_>0) && (failed_==false) );
00392 
00393   ordered.resize( nodes_ );
00394   int loc = 0;
00395   for( int i = 0; i < depth_; ++i )
00396   {
00397     int currLevel = depth_ - (i+1);
00398     int currWidth = levelSets_[currLevel].size();
00399     for( int j = 0; j < currWidth; ++j )
00400       ordered[loc++] = levelSets_[currLevel][currWidth-j-1];
00401   }
00402 }
00403 
00404 } //namespace EpetraExt
00405 #endif //HAVE_GRAPH_REORDERINGS
00406 #endif //HAVE_EXPERIMENTAL

Generated on Wed May 12 21:24:46 2010 for EpetraExt by  doxygen 1.4.7