EpetraExt_MapColoring.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_MapColoring.h>
00030 
00031 #include <EpetraExt_Transpose_CrsGraph.h>
00032 #include <EpetraExt_Overlap_CrsGraph.h>
00033 
00034 #include <Epetra_CrsGraph.h>
00035 #include <Epetra_IntVector.h>
00036 #include <Epetra_MapColoring.h>
00037 #include <Epetra_Map.h>
00038 #include <Epetra_Comm.h>
00039 #include <Epetra_Util.h>
00040 #include <Epetra_Import.h>
00041 #include <Epetra_Export.h>
00042 
00043 #include <Epetra_Time.h>
00044 
00045 #include <vector>
00046 #include <set>
00047 #include <map>
00048 
00049 using std::vector;
00050 using std::set;
00051 using std::map;
00052 
00053 namespace EpetraExt {
00054 
00055 CrsGraph_MapColoring::NewTypeRef
00056 CrsGraph_MapColoring::
00057 operator()( OriginalTypeRef orig  )
00058 {
00059   Epetra_Time timer( orig.Comm() );
00060 
00061   origObj_ = &orig;
00062 
00063   const Epetra_BlockMap & RowMap = orig.RowMap();
00064   int nRows = RowMap.NumMyElements();
00065   const Epetra_BlockMap & ColMap = orig.ColMap();
00066   int nCols = ColMap.NumMyElements();
00067 
00068   int MyPID = RowMap.Comm().MyPID();
00069 
00070   if( verbosity_ > 1 ) cout << "RowMap:\n" << RowMap;
00071   if( verbosity_ > 1 ) cout << "ColMap:\n" << ColMap;
00072 
00073   Epetra_MapColoring * ColorMap = 0;
00074 
00075   if( algo_ == GREEDY || algo_ == LUBY )
00076   {
00077 
00078     Epetra_CrsGraph * base = &( const_cast<Epetra_CrsGraph&>(orig) );
00079 
00080     int NumIndices;
00081     int * Indices;
00082 
00083     double wTime1 = timer.WallTime();
00084 
00085     // For parallel applications, add in boundaries to coloring
00086     bool distributedGraph = RowMap.DistributedGlobal();
00087     if( distributedGraph )
00088     {
00089       base = new Epetra_CrsGraph( Copy, ColMap, ColMap, 0 );
00090 
00091       for( int i = 0; i < nRows; ++i )
00092       {
00093         orig.ExtractMyRowView( i, NumIndices, Indices );
00094         base->InsertMyIndices( i, NumIndices, Indices );
00095 
00096         //Do this with a single insert
00097         //Is this the right thing?
00098         for( int j = 0; j < NumIndices; ++j )
00099           if( Indices[j] >= nRows )
00100             base->InsertMyIndices( Indices[j], 1, &i );
00101       } 
00102       base->FillComplete();
00103     }
00104 
00105     if( verbosity_ > 1 ) cout << "Base Graph:\n" << *base << endl;
00106 
00107     double wTime2 = timer.WallTime();
00108     if( verbosity_ > 0 )
00109       cout << "EpetraExt::MapColoring [INSERT BOUNDARIES] Time: " << wTime2-wTime1 << endl;
00110 
00111     //Generate Local Distance-1 Adjacency Graph
00112     //(Transpose of orig excluding non-local column indices)
00113     EpetraExt::CrsGraph_Transpose transposeTransform( true );
00114     Epetra_CrsGraph & Adj1 = transposeTransform( *base );
00115     if( verbosity_ > 1 ) cout << "Adjacency 1 Graph!\n" << Adj1;
00116 
00117     wTime1 = timer.WallTime();
00118     if( verbosity_ > 0 )
00119       cout << "EpetraExt::MapColoring [TRANSPOSE GRAPH] Time: " << wTime1-wTime2 << endl;
00120 
00121     int Delta = Adj1.MaxNumIndices();
00122     if( verbosity_ > 0 ) cout << endl << "Delta: " << Delta << endl;
00123 
00124     //Generation of Local Distance-2 Adjacency Graph
00125     Epetra_CrsGraph * Adj2 = &Adj1;
00126     if( !distance1_ )
00127     {
00128       Adj2 = new Epetra_CrsGraph( Copy, ColMap, ColMap, 0 );
00129       int NumAdj1Indices;
00130       int * Adj1Indices;
00131       for( int i = 0; i < nCols; ++i )
00132       {
00133         Adj1.ExtractMyRowView( i, NumAdj1Indices, Adj1Indices );
00134 
00135         set<int> Cols;
00136         for( int j = 0; j < NumAdj1Indices; ++j )
00137         {
00138           base->ExtractMyRowView( Adj1Indices[j], NumIndices, Indices );
00139           for( int k = 0; k < NumIndices; ++k )
00140             if( Indices[k] < nCols ) Cols.insert( Indices[k] );
00141         }
00142         int nCols2 = Cols.size();
00143         vector<int> ColVec( nCols2 );
00144         set<int>::iterator iterIS = Cols.begin();
00145         set<int>::iterator iendIS = Cols.end();
00146         for( int j = 0 ; iterIS != iendIS; ++iterIS, ++j ) ColVec[j] = *iterIS;
00147         Adj2->InsertMyIndices( i, nCols2, &ColVec[0] );
00148       }
00149       Adj2->FillComplete();
00150 
00151       if( verbosity_ > 1 ) cout << "Adjacency 2 Graph!\n" << *Adj2;
00152 
00153       wTime2 = timer.WallTime();
00154       if( verbosity_ > 0 )
00155         cout << "EpetraExt::MapColoring [GEN DIST-2 GRAPH] Time: " << wTime2-wTime1 << endl;
00156     }
00157 
00158     wTime2 = timer.WallTime();
00159 
00160     ColorMap = new Epetra_MapColoring( ColMap );
00161 
00162     vector<int> rowOrder( nCols );
00163     if( reordering_ == 0 || reordering_ == 1 ) 
00164     {
00165       multimap<int,int> adjMap;
00166       typedef multimap<int,int>::value_type adjMapValueType;
00167       for( int i = 0; i < nCols; ++i )
00168         adjMap.insert( adjMapValueType( Adj2->NumMyIndices(i), i ) );
00169       multimap<int,int>::iterator iter = adjMap.begin();
00170       multimap<int,int>::iterator end = adjMap.end();
00171       if( reordering_ == 0 ) //largest first (less colors)
00172       {
00173         for( int i = 1; iter != end; ++iter, ++i )
00174           rowOrder[ nCols - i ] = (*iter).second;
00175       }
00176       else                  //smallest first (better balance)
00177       {
00178         for( int i = 0; iter != end; ++iter, ++i )
00179           rowOrder[ i ] = (*iter).second;
00180       }
00181     }
00182     else if( reordering_ == 2 ) //random
00183     {
00184       for( int i = 0; i < nCols; ++i )
00185         rowOrder[ i ] = i;
00186 #ifndef TFLOP
00187       random_shuffle( rowOrder.begin(), rowOrder.end() );
00188 #endif
00189     }
00190 
00191     wTime1 = timer.WallTime();
00192     if( verbosity_ > 0 )
00193       cout << "EpetraExt::MapColoring [REORDERING] Time: " << wTime1-wTime2 << endl;
00194 
00195     //Application of Greedy Algorithm to generate Color Map
00196     if( algo_ == GREEDY )
00197     {
00198       for( int col = 0; col < nCols; ++col )
00199       {
00200         Adj2->ExtractMyRowView( rowOrder[col], NumIndices, Indices );
00201 
00202         set<int> usedColors;
00203         int color;
00204         for( int i = 0; i < NumIndices; ++i )
00205         {
00206           color = (*ColorMap)[ Indices[i] ];
00207           if( color > 0 ) usedColors.insert( color );
00208           color = 0;
00209           int testcolor = 1;
00210           while( !color )
00211           {
00212             if( !usedColors.count( testcolor ) ) color = testcolor;
00213             ++testcolor;
00214           }
00215         }
00216         (*ColorMap)[ rowOrder[col] ] = color;
00217       }
00218 
00219       wTime2 = timer.WallTime();
00220       if( verbosity_ > 0 )
00221         cout << "EpetraExt::MapColoring [GREEDY COLORING] Time: " << wTime2-wTime1 << endl;
00222       if( verbosity_ > 0 )
00223         cout << "Num GREEDY Colors: " << ColorMap->NumColors() << endl;
00224     }
00225     else if( algo_ == LUBY )
00226     {
00227        //Assign Random Keys To Rows
00228       Epetra_Util util;
00229       vector<int> Keys(nCols);
00230       vector<int> State(nCols,-1);
00231 
00232       for( int col = 0; col < nCols; ++col )
00233         Keys[col] = util.RandomInt();
00234 
00235       int NumRemaining = nCols;
00236       int CurrentColor = 1;
00237 
00238       while( NumRemaining > 0 )
00239       {
00240         //maximal independent set
00241         while( NumRemaining > 0 )
00242         {
00243           NumRemaining = 0;
00244 
00245           //zero out everyone less than neighbor
00246           for( int i = 0; i < nCols; ++i )
00247           {
00248             int col = rowOrder[i];
00249             if( State[col] < 0 )
00250             {
00251               Adj2->ExtractMyRowView( col, NumIndices, Indices );
00252               int MyKey = Keys[col];
00253               for( int j = 0; j < NumIndices; ++j )
00254                 if( col != Indices[j] && State[ Indices[j] ] < 0 )
00255                 {
00256                   if( MyKey > Keys[ Indices[j] ] ) State[ Indices[j] ] = 0;
00257                   else                             State[ col ] = 0;
00258                 }
00259             }
00260           }
00261 
00262           //assign -1's to current color
00263           for( int col = 0; col < nCols; ++col )
00264           {
00265             if( State[col] < 0 )
00266               State[col] = CurrentColor;
00267           }
00268 
00269           //reinstate any zero not neighboring current color
00270           for( int col = 0; col < nCols; ++col )
00271             if( State[col] == 0 )
00272             {
00273               Adj2->ExtractMyRowView( col, NumIndices, Indices );
00274               bool flag = false;
00275               for( int i = 0; i < NumIndices; ++i )
00276                 if( col != Indices[i] && State[ Indices[i] ] == CurrentColor )
00277                 {
00278                   flag = true;
00279                   break;
00280                 }
00281               if( !flag )
00282               {
00283                 State[col] = -1;
00284                 ++NumRemaining;
00285               }
00286             }
00287         }
00288 
00289         //Reset Status for all non-colored nodes
00290         for( int col = 0; col < nCols; ++col )
00291           if( State[col] == 0 )
00292           {
00293             State[col] = -1;
00294             ++NumRemaining;
00295           }
00296 
00297         if( verbosity_ > 2 )
00298         {
00299           cout << "Finished Color: " << CurrentColor << endl;
00300           cout << "NumRemaining: " << NumRemaining << endl;
00301         }
00302 
00303         //New color
00304         ++CurrentColor;
00305       }
00306 
00307       for( int col = 0; col < nCols; ++col )
00308         (*ColorMap)[col] = State[col]-1;
00309 
00310       wTime2 = timer.WallTime();
00311       if( verbosity_ > 0 )
00312         cout << "EpetraExt::MapColoring [LUBI COLORING] Time: " << wTime2-wTime1 << endl;
00313       if( verbosity_ > 0 )
00314         cout << "Num LUBI Colors: " << ColorMap->NumColors() << endl;
00315 
00316     }
00317     else
00318       abort(); //UNKNOWN ALGORITHM
00319 
00320     if( distributedGraph ) delete base;
00321     if( !distance1_ ) delete Adj2;
00322   }
00323   else
00324   {
00325     //Generate Parallel Adjacency-2 Graph
00326     EpetraExt::CrsGraph_Overlap OverlapTrans(1);
00327     Epetra_CrsGraph & OverlapGraph = OverlapTrans( orig );
00328 
00329     if( verbosity_ > 1 ) cout << "OverlapGraph:\n" << OverlapGraph;
00330 
00331     Epetra_CrsGraph * Adj2 = &orig;
00332 
00333     int NumIndices;
00334     int * Indices;
00335     if( !distance1_ )
00336     {
00337       Adj2 = new Epetra_CrsGraph( Copy, RowMap, OverlapGraph.ColMap(), 0 );
00338       int NumAdj1Indices;
00339       int * Adj1Indices;
00340       for( int i = 0; i < nRows; ++i )
00341       {
00342         OverlapGraph.ExtractMyRowView( i, NumAdj1Indices, Adj1Indices );
00343 
00344         set<int> Cols;
00345         for( int j = 0; j < NumAdj1Indices; ++j )
00346         {
00347           int GID = OverlapGraph.LRID( OverlapGraph.GCID( Adj1Indices[j] ) );
00348           OverlapGraph.ExtractMyRowView( GID, NumIndices, Indices );
00349           for( int k = 0; k < NumIndices; ++k ) Cols.insert( Indices[k] );
00350         }
00351         int nCols2 = Cols.size();
00352         vector<int> ColVec( nCols2 );
00353         set<int>::iterator iterIS = Cols.begin();
00354         set<int>::iterator iendIS = Cols.end();
00355         for( int j = 0 ; iterIS != iendIS; ++iterIS, ++j ) ColVec[j] = *iterIS;
00356         Adj2->InsertMyIndices( i, nCols2, &ColVec[0] );
00357       }
00358       int flag = Adj2->FillComplete();
00359       assert( flag == 0 );
00360       RowMap.Comm().Barrier();
00361       if( verbosity_ > 1 ) cout << "Adjacency 2 Graph!\n" << *Adj2;
00362     }
00363 
00364     //collect GIDs on boundary
00365     vector<int> boundaryGIDs;
00366     vector<int> interiorGIDs;
00367     for( int row = 0; row < nRows; ++row )
00368     {
00369       Adj2->ExtractMyRowView( row, NumIndices, Indices );
00370       bool testFlag = false;
00371       for( int i = 0; i < NumIndices; ++i )
00372         if( Indices[i] >= nRows ) testFlag = true;
00373       if( testFlag ) boundaryGIDs.push_back( Adj2->GRID(row) );
00374       else           interiorGIDs.push_back( Adj2->GRID(row) );
00375     }
00376 
00377     int LocalBoundarySize = boundaryGIDs.size();
00378 
00379     Epetra_Map BoundaryMap( -1, boundaryGIDs.size(), &boundaryGIDs[0], 0, RowMap.Comm() );
00380     if( verbosity_ > 1 ) cout << "BoundaryMap:\n" << BoundaryMap;
00381     
00382     int BoundarySize = BoundaryMap.NumGlobalElements();
00383     Epetra_MapColoring BoundaryColoring( BoundaryMap );
00384 
00385     if( algo_ == PSEUDO_PARALLEL )
00386     {
00387       Epetra_Map BoundaryIndexMap( BoundarySize, LocalBoundarySize, 0, RowMap.Comm() );
00388       if( verbosity_ > 1) cout << "BoundaryIndexMap:\n" << BoundaryIndexMap;
00389 
00390       Epetra_IntVector bGIDs( View, BoundaryIndexMap, &boundaryGIDs[0] );
00391       if( verbosity_ > 1) cout << "BoundaryGIDs:\n" << bGIDs;
00392 
00393       int NumLocalBs = 0;
00394       if( !RowMap.Comm().MyPID() ) NumLocalBs = BoundarySize;
00395      
00396       Epetra_Map LocalBoundaryIndexMap( BoundarySize, NumLocalBs, 0, RowMap.Comm() );
00397       if( verbosity_ > 1) cout << "LocalBoundaryIndexMap:\n" << LocalBoundaryIndexMap;
00398 
00399       Epetra_IntVector lbGIDs( LocalBoundaryIndexMap );
00400       Epetra_Import lbImport( LocalBoundaryIndexMap, BoundaryIndexMap );
00401       lbGIDs.Import( bGIDs, lbImport, Insert );
00402       if( verbosity_ > 1) cout << "LocalBoundaryGIDs:\n" << lbGIDs;
00403 
00404       Epetra_Map LocalBoundaryMap( BoundarySize, NumLocalBs, lbGIDs.Values(), 0, RowMap.Comm() );
00405       if( verbosity_ > 1) cout << "LocalBoundaryMap:\n" << LocalBoundaryMap;
00406 
00407       Epetra_CrsGraph LocalBoundaryGraph( Copy, LocalBoundaryMap, LocalBoundaryMap, 0 );
00408       Epetra_Import LocalBoundaryImport( LocalBoundaryMap, Adj2->RowMap() );
00409       LocalBoundaryGraph.Import( *Adj2, LocalBoundaryImport, Insert );
00410       LocalBoundaryGraph.FillComplete();
00411       if( verbosity_ > 1 ) cout << "LocalBoundaryGraph:\n " << LocalBoundaryGraph;
00412 
00413       EpetraExt::CrsGraph_MapColoring BoundaryTrans( GREEDY, reordering_, distance1_, verbosity_ );
00414       Epetra_MapColoring & LocalBoundaryColoring = BoundaryTrans( LocalBoundaryGraph );
00415       if( verbosity_ > 1 ) cout << "LocalBoundaryColoring:\n " << LocalBoundaryColoring;
00416 
00417       Epetra_Export BoundaryExport( LocalBoundaryMap, BoundaryMap );
00418       BoundaryColoring.Export( LocalBoundaryColoring, BoundaryExport, Insert );
00419     }
00420     else if( algo_ == JONES_PLASSMAN )
00421     {
00422     /* Alternative Distrib. Memory Coloring of Boundary based on JonesPlassman(sic) paper
00423      * 1.Random number assignment to all boundary nodes using GID as seed to function
00424      * (This allows any processor to compute adj. off proc values with a local computation)
00425      * 2.Find all nodes greater than any neighbor off processor, color them.
00426      * 3.Send colored node info to neighbors
00427      * 4.Constrained color all nodes with all off proc neighbors smaller or colored.
00428      * 5.Goto 3
00429      */
00430 
00431       vector<int> OverlapBoundaryGIDs( boundaryGIDs );
00432       for( int i = nRows; i < Adj2->ColMap().NumMyElements(); ++i )
00433         OverlapBoundaryGIDs.push_back( Adj2->ColMap().GID(i) );
00434 
00435       int OverlapBoundarySize = OverlapBoundaryGIDs.size();
00436       Epetra_Map BoundaryColMap( -1, OverlapBoundarySize, &OverlapBoundaryGIDs[0], 0, RowMap.Comm() );
00437 
00438       Epetra_CrsGraph BoundaryGraph( Copy, BoundaryMap, BoundaryColMap, 0 );
00439       Epetra_Import BoundaryImport( BoundaryMap, Adj2->RowMap() );
00440       BoundaryGraph.Import( *Adj2, BoundaryImport, Insert );
00441       BoundaryGraph.FillComplete();
00442       if( verbosity_ > 1) cout << "BoundaryGraph:\n" << BoundaryGraph;
00443 
00444       Epetra_Import ReverseOverlapBoundaryImport( BoundaryMap, BoundaryColMap );
00445       Epetra_Import OverlapBoundaryImport( BoundaryColMap, BoundaryMap );
00446 
00447       int Colored = 0;
00448       int GlobalColored = 0;
00449       int Level = 0;
00450       Epetra_MapColoring OverlapBoundaryColoring( BoundaryColMap );
00451 
00452       //Setup random integers for boundary nodes
00453       Epetra_IntVector BoundaryValues( BoundaryMap );
00454       Epetra_Util Util;
00455       Util.SetSeed( 47954118 * (MyPID+1) );
00456       for( int i=0; i < LocalBoundarySize; ++i )
00457       {
00458         int val = Util.RandomInt();
00459         if( val < 0 ) val *= -1;
00460         BoundaryValues[i] = val;
00461       }
00462 
00463       //Get Random Values for External Boundary
00464       Epetra_IntVector OverlapBoundaryValues( BoundaryColMap );
00465       OverlapBoundaryValues.Import( BoundaryValues, OverlapBoundaryImport, Insert );
00466 
00467       while( GlobalColored < BoundarySize )
00468       {
00469   //Find current "Level" of boundary indices to color
00470   int NumIndices;
00471   int * Indices;
00472   vector<int> LevelIndices;
00473   for( int i = 0; i < LocalBoundarySize; ++i )
00474   {
00475           if( !OverlapBoundaryColoring[i] )
00476           {
00477             //int MyVal = PRAND(BoundaryColMap.GID(i));
00478             int MyVal = OverlapBoundaryValues[i];
00479             BoundaryGraph.ExtractMyRowView( i, NumIndices, Indices );
00480       bool ColorFlag = true;
00481       int Loc = 0;
00482       while( Loc<NumIndices && Indices[Loc]<LocalBoundarySize ) ++Loc;
00483       for( int j = Loc; j < NumIndices; ++j )
00484               if( (OverlapBoundaryValues[Indices[j]]>MyVal)
00485             && !OverlapBoundaryColoring[Indices[j]] )
00486               {
00487                 ColorFlag = false;
00488     break;
00489               }
00490             if( ColorFlag ) LevelIndices.push_back(i);
00491           }
00492         }
00493 
00494   if( verbosity_ > 1 )
00495         {
00496           cout << MyPID << " Level Indices: ";
00497     int Lsize = (int) LevelIndices.size();
00498     for( int i = 0; i < Lsize; ++i ) cout << LevelIndices[i] << " ";
00499     cout << endl;
00500         }
00501 
00502         //Greedy coloring of current level
00503   set<int> levelColors;
00504   int Lsize = (int) LevelIndices.size();
00505         for( int i = 0; i < Lsize; ++i )
00506         {
00507           BoundaryGraph.ExtractMyRowView( LevelIndices[i], NumIndices, Indices );
00508 
00509           set<int> usedColors;
00510           int color;
00511           for( int j = 0; j < NumIndices; ++j )
00512           {
00513             color = OverlapBoundaryColoring[ Indices[j] ];
00514             if( color > 0 ) usedColors.insert( color );
00515             color = 0;
00516             int testcolor = 1;
00517             while( !color )
00518             {
00519               if( !usedColors.count( testcolor ) ) color = testcolor;
00520               ++testcolor;
00521             }
00522           }
00523           OverlapBoundaryColoring[ LevelIndices[i] ] = color;
00524     levelColors.insert( color );
00525         }
00526 
00527   if( verbosity_ > 2 ) cout << MyPID << " Level: " << Level << " Count: " << LevelIndices.size() << " NumColors: " << levelColors.size() << endl;
00528 
00529   if( verbosity_ > 2 ) cout << "Current Level Boundary Coloring:\n" << OverlapBoundaryColoring;
00530 
00531   //Update off processor coloring info
00532   BoundaryColoring.Import( OverlapBoundaryColoring, ReverseOverlapBoundaryImport, Insert );
00533   OverlapBoundaryColoring.Import( BoundaryColoring, OverlapBoundaryImport, Insert );
00534         Colored += LevelIndices.size();
00535   Level++;
00536 
00537   RowMap.Comm().SumAll( &Colored, &GlobalColored, 1 );
00538   if( verbosity_ > 2 ) cout << "Num Globally Colored: " << GlobalColored << " from Num Global Boundary Nodes: " << BoundarySize << endl;
00539       }
00540     }
00541 
00542     if( verbosity_ > 1 ) cout << "BoundaryColoring:\n " << BoundaryColoring;
00543 
00544     Epetra_MapColoring RowColorMap( RowMap );
00545 
00546     //Add Boundary Colors
00547     for( int i = 0; i < LocalBoundarySize; ++i )
00548     {
00549       int GID = BoundaryMap.GID(i);
00550       RowColorMap(GID) = BoundaryColoring(GID);
00551     }
00552 
00553     Epetra_MapColoring Adj2ColColorMap( Adj2->ColMap() );
00554     Epetra_Import Adj2Import( Adj2->ColMap(), RowMap );
00555     Adj2ColColorMap.Import( RowColorMap, Adj2Import, Insert );
00556 
00557     if( verbosity_ > 1 ) cout << "RowColoringMap:\n " << RowColorMap;
00558     if( verbosity_ > 1 ) cout << "Adj2ColColorMap:\n " << Adj2ColColorMap;
00559 
00560     vector<int> rowOrder( nRows );
00561     if( reordering_ == 0 || reordering_ == 1 ) 
00562     {
00563       multimap<int,int> adjMap;
00564       typedef multimap<int,int>::value_type adjMapValueType;
00565       for( int i = 0; i < nRows; ++i )
00566         adjMap.insert( adjMapValueType( Adj2->NumMyIndices(i), i ) );
00567       multimap<int,int>::iterator iter = adjMap.begin();
00568       multimap<int,int>::iterator end = adjMap.end();
00569       if( reordering_ == 0 ) //largest first (less colors)
00570       {
00571         for( int i = 1; iter != end; ++iter, ++i )
00572           rowOrder[nRows-i] = (*iter).second;
00573       }
00574       else                  //smallest first (better balance)
00575       {
00576         for( int i = 0; iter != end; ++iter, ++i )
00577           rowOrder[i] = (*iter).second;
00578       }
00579     }
00580     else if( reordering_ == 2 ) //random
00581     {
00582       for( int i = 0; i < nRows; ++i )
00583         rowOrder[i] = i;
00584 #ifdef TFLOP
00585       random_shuffle( rowOrder.begin(), rowOrder.end() );
00586 #endif
00587     }
00588 
00589     //Constrained greedy coloring of interior
00590     set<int> InteriorColors;
00591     for( int row = 0; row < nRows; ++row )
00592     {
00593       if( !RowColorMap[ rowOrder[row] ] )
00594       {
00595         Adj2->ExtractMyRowView( rowOrder[row], NumIndices, Indices );
00596 
00597         set<int> usedColors;
00598         int color;
00599         for( int i = 0; i < NumIndices; ++i )
00600         {
00601           color = Adj2ColColorMap[ Indices[i] ];
00602           if( color > 0 ) usedColors.insert( color );
00603           color = 0;
00604           int testcolor = 1;
00605           while( !color )
00606           {
00607             if( !usedColors.count( testcolor ) ) color = testcolor;
00608             ++testcolor;
00609           }
00610         }
00611         Adj2ColColorMap[ rowOrder[row] ] = color;
00612   InteriorColors.insert( color );
00613       }
00614     }
00615     if( verbosity_ > 1 ) cout << MyPID << " Num Interior Colors: " << InteriorColors.size() << endl;
00616     if( verbosity_ > 1 ) cout << "RowColorMap after Greedy:\n " << RowColorMap;
00617 
00618     ColorMap = new Epetra_MapColoring( ColMap );
00619     Epetra_Import ColImport( ColMap, Adj2->ColMap() );
00620     ColorMap->Import( Adj2ColColorMap, ColImport, Insert );
00621 
00622     if( !distance1_ ) delete Adj2;
00623   }
00624 
00625   if( verbosity_ > 0 ) cout << MyPID << " ColorMap Color Count: " << ColorMap->NumColors() << endl;
00626   if( verbosity_ > 1 ) cout << "ColorMap!\n" << *ColorMap;
00627 
00628   newObj_ = ColorMap;
00629 
00630   return *ColorMap;
00631 }
00632 
00633 } // namespace EpetraExt

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