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

Generated on Thu Sep 18 12:31:44 2008 for EpetraExt by doxygen 1.3.9.1