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