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