EpetraExt Package Browser (Single Doxygen Collection) 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 (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 ) std::cout << "RowMap:\n" << RowMap;
00071   if( verbosity_ > 1 ) std::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 ) std::cout << "Base Graph:\n" << *base << std::endl;
00106 
00107     double wTime2 = timer.WallTime();
00108     if( verbosity_ > 0 )
00109       std::cout << "EpetraExt::MapColoring [INSERT BOUNDARIES] Time: " << wTime2-wTime1 << std::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 ) std::cout << "Adjacency 1 Graph!\n" << Adj1;
00116 
00117     wTime1 = timer.WallTime();
00118     if( verbosity_ > 0 )
00119       std::cout << "EpetraExt::MapColoring [TRANSPOSE GRAPH] Time: " << wTime1-wTime2 << std::endl;
00120 
00121     int Delta = Adj1.MaxNumIndices();
00122     if( verbosity_ > 0 ) std::cout << std::endl << "Delta: " << Delta << std::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         std::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 ) std::cout << "Adjacency 2 Graph!\n" << *Adj2;
00152 
00153       wTime2 = timer.WallTime();
00154       if( verbosity_ > 0 )
00155         std::cout << "EpetraExt::MapColoring [GEN DIST-2 GRAPH] Time: " << wTime2-wTime1 << std::endl;
00156     }
00157 
00158     wTime2 = timer.WallTime();
00159 
00160     ColorMap = new Epetra_MapColoring( ColMap );
00161 
00162     std::vector<int> rowOrder( nCols );
00163     if( reordering_ == 0 || reordering_ == 1 ) 
00164     {
00165       std::multimap<int,int> adjMap;
00166       typedef std::multimap<int,int>::value_type adjMapValueType;
00167       for( int i = 0; i < nCols; ++i )
00168         adjMap.insert( adjMapValueType( Adj2->NumMyIndices(i), i ) );
00169       std::multimap<int,int>::iterator iter = adjMap.begin();
00170       std::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       std::random_shuffle( rowOrder.begin(), rowOrder.end() );
00188 #endif
00189     }
00190 
00191     wTime1 = timer.WallTime();
00192     if( verbosity_ > 0 )
00193       std::cout << "EpetraExt::MapColoring [REORDERING] Time: " << wTime1-wTime2 << std::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         std::cout << "EpetraExt::MapColoring [GREEDY COLORING] Time: " << wTime2-wTime1 << std::endl;
00222       if( verbosity_ > 0 )
00223         std::cout << "Num GREEDY Colors: " << ColorMap->NumColors() << std::endl;
00224     }
00225     else if( algo_ == LUBY )
00226     {
00227        //Assign Random Keys To Rows
00228       Epetra_Util util;
00229       std::vector<int> Keys(nCols);
00230       std::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           std::cout << "Finished Color: " << CurrentColor << std::endl;
00300           std::cout << "NumRemaining: " << NumRemaining << std::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         std::cout << "EpetraExt::MapColoring [LUBI COLORING] Time: " << wTime2-wTime1 << std::endl;
00313       if( verbosity_ > 0 )
00314         std::cout << "Num LUBI Colors: " << ColorMap->NumColors() << std::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 ) std::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         std::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 ) std::cout << "Adjacency 2 Graph!\n" << *Adj2;
00362     }
00363 
00364     //collect GIDs on boundary
00365     std::vector<int> boundaryGIDs;
00366     std::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(),
00380       LocalBoundarySize ? &boundaryGIDs[0]: 0,
00381       0, RowMap.Comm() );
00382     if( verbosity_ > 1 ) std::cout << "BoundaryMap:\n" << BoundaryMap;
00383     
00384     int BoundarySize = BoundaryMap.NumGlobalElements();
00385     Epetra_MapColoring BoundaryColoring( BoundaryMap );
00386 
00387     if( algo_ == PSEUDO_PARALLEL )
00388     {
00389       Epetra_Map BoundaryIndexMap( BoundarySize, LocalBoundarySize, 0, RowMap.Comm() );
00390       if( verbosity_ > 1) std::cout << "BoundaryIndexMap:\n" << BoundaryIndexMap;
00391 
00392       Epetra_IntVector bGIDs( View, BoundaryIndexMap, &boundaryGIDs[0] );
00393       if( verbosity_ > 1) std::cout << "BoundaryGIDs:\n" << bGIDs;
00394 
00395       int NumLocalBs = 0;
00396       if( !RowMap.Comm().MyPID() ) NumLocalBs = BoundarySize;
00397      
00398       Epetra_Map LocalBoundaryIndexMap( BoundarySize, NumLocalBs, 0, RowMap.Comm() );
00399       if( verbosity_ > 1) std::cout << "LocalBoundaryIndexMap:\n" << LocalBoundaryIndexMap;
00400 
00401       Epetra_IntVector lbGIDs( LocalBoundaryIndexMap );
00402       Epetra_Import lbImport( LocalBoundaryIndexMap, BoundaryIndexMap );
00403       lbGIDs.Import( bGIDs, lbImport, Insert );
00404       if( verbosity_ > 1) std::cout << "LocalBoundaryGIDs:\n" << lbGIDs;
00405 
00406       Epetra_Map LocalBoundaryMap( BoundarySize, NumLocalBs, lbGIDs.Values(), 0, RowMap.Comm() );
00407       if( verbosity_ > 1) std::cout << "LocalBoundaryMap:\n" << LocalBoundaryMap;
00408 
00409       Epetra_CrsGraph LocalBoundaryGraph( Copy, LocalBoundaryMap, LocalBoundaryMap, 0 );
00410       Epetra_Import LocalBoundaryImport( LocalBoundaryMap, Adj2->RowMap() );
00411       LocalBoundaryGraph.Import( *Adj2, LocalBoundaryImport, Insert );
00412       LocalBoundaryGraph.FillComplete();
00413       if( verbosity_ > 1 ) std::cout << "LocalBoundaryGraph:\n " << LocalBoundaryGraph;
00414 
00415       EpetraExt::CrsGraph_MapColoring BoundaryTrans( GREEDY, reordering_, distance1_, verbosity_ );
00416       Epetra_MapColoring & LocalBoundaryColoring = BoundaryTrans( LocalBoundaryGraph );
00417       if( verbosity_ > 1 ) std::cout << "LocalBoundaryColoring:\n " << LocalBoundaryColoring;
00418 
00419       Epetra_Export BoundaryExport( LocalBoundaryMap, BoundaryMap );
00420       BoundaryColoring.Export( LocalBoundaryColoring, BoundaryExport, Insert );
00421     }
00422     else if( algo_ == JONES_PLASSMAN )
00423     {
00424     /* Alternative Distrib. Memory Coloring of Boundary based on JonesPlassman(sic) paper
00425      * 1.Random number assignment to all boundary nodes using GID as seed to function
00426      * (This allows any processor to compute adj. off proc values with a local computation)
00427      * 2.Find all nodes greater than any neighbor off processor, color them.
00428      * 3.Send colored node info to neighbors
00429      * 4.Constrained color all nodes with all off proc neighbors smaller or colored.
00430      * 5.Goto 3
00431      */
00432 
00433       std::vector<int> OverlapBoundaryGIDs( boundaryGIDs );
00434       for( int i = nRows; i < Adj2->ColMap().NumMyElements(); ++i )
00435         OverlapBoundaryGIDs.push_back( Adj2->ColMap().GID(i) );
00436 
00437       int OverlapBoundarySize = OverlapBoundaryGIDs.size();
00438       Epetra_Map BoundaryColMap( -1, OverlapBoundarySize,
00439         OverlapBoundarySize ? &OverlapBoundaryGIDs[0] : 0,
00440         0, RowMap.Comm() );
00441 
00442       Epetra_CrsGraph BoundaryGraph( Copy, BoundaryMap, BoundaryColMap, 0 );
00443       Epetra_Import BoundaryImport( BoundaryMap, Adj2->RowMap() );
00444       BoundaryGraph.Import( *Adj2, BoundaryImport, Insert );
00445       BoundaryGraph.FillComplete();
00446       if( verbosity_ > 1) std::cout << "BoundaryGraph:\n" << BoundaryGraph;
00447 
00448       Epetra_Import ReverseOverlapBoundaryImport( BoundaryMap, BoundaryColMap );
00449       Epetra_Import OverlapBoundaryImport( BoundaryColMap, BoundaryMap );
00450 
00451       int Colored = 0;
00452       int GlobalColored = 0;
00453       int Level = 0;
00454       Epetra_MapColoring OverlapBoundaryColoring( BoundaryColMap );
00455 
00456       //Setup random integers for boundary nodes
00457       Epetra_IntVector BoundaryValues( BoundaryMap );
00458       Epetra_Util Util;
00459       Util.SetSeed( 47954118 * (MyPID+1) );
00460       for( int i=0; i < LocalBoundarySize; ++i )
00461       {
00462         int val = Util.RandomInt();
00463         if( val < 0 ) val *= -1;
00464         BoundaryValues[i] = val;
00465       }
00466 
00467       //Get Random Values for External Boundary
00468       Epetra_IntVector OverlapBoundaryValues( BoundaryColMap );
00469       OverlapBoundaryValues.Import( BoundaryValues, OverlapBoundaryImport, Insert );
00470 
00471       while( GlobalColored < BoundarySize )
00472       {
00473   //Find current "Level" of boundary indices to color
00474   int NumIndices;
00475   int * Indices;
00476   std::vector<int> LevelIndices;
00477   for( int i = 0; i < LocalBoundarySize; ++i )
00478   {
00479           if( !OverlapBoundaryColoring[i] )
00480           {
00481             //int MyVal = PRAND(BoundaryColMap.GID(i));
00482             int MyVal = OverlapBoundaryValues[i];
00483             BoundaryGraph.ExtractMyRowView( i, NumIndices, Indices );
00484       bool ColorFlag = true;
00485       int Loc = 0;
00486       while( Loc<NumIndices && Indices[Loc]<LocalBoundarySize ) ++Loc;
00487       for( int j = Loc; j < NumIndices; ++j )
00488               if( (OverlapBoundaryValues[Indices[j]]>MyVal)
00489             && !OverlapBoundaryColoring[Indices[j]] )
00490               {
00491                 ColorFlag = false;
00492     break;
00493               }
00494             if( ColorFlag ) LevelIndices.push_back(i);
00495           }
00496         }
00497 
00498   if( verbosity_ > 1 )
00499         {
00500           std::cout << MyPID << " Level Indices: ";
00501     int Lsize = (int) LevelIndices.size();
00502     for( int i = 0; i < Lsize; ++i ) std::cout << LevelIndices[i] << " ";
00503     std::cout << std::endl;
00504         }
00505 
00506         //Greedy coloring of current level
00507   set<int> levelColors;
00508   int Lsize = (int) LevelIndices.size();
00509         for( int i = 0; i < Lsize; ++i )
00510         {
00511           BoundaryGraph.ExtractMyRowView( LevelIndices[i], NumIndices, Indices );
00512 
00513           set<int> usedColors;
00514           int color;
00515           for( int j = 0; j < NumIndices; ++j )
00516           {
00517             color = OverlapBoundaryColoring[ Indices[j] ];
00518             if( color > 0 ) usedColors.insert( color );
00519             color = 0;
00520             int testcolor = 1;
00521             while( !color )
00522             {
00523               if( !usedColors.count( testcolor ) ) color = testcolor;
00524               ++testcolor;
00525             }
00526           }
00527           OverlapBoundaryColoring[ LevelIndices[i] ] = color;
00528     levelColors.insert( color );
00529         }
00530 
00531   if( verbosity_ > 2 ) std::cout << MyPID << " Level: " << Level << " Count: " << LevelIndices.size() << " NumColors: " << levelColors.size() << std::endl;
00532 
00533   if( verbosity_ > 2 ) std::cout << "Current Level Boundary Coloring:\n" << OverlapBoundaryColoring;
00534 
00535   //Update off processor coloring info
00536   BoundaryColoring.Import( OverlapBoundaryColoring, ReverseOverlapBoundaryImport, Insert );
00537   OverlapBoundaryColoring.Import( BoundaryColoring, OverlapBoundaryImport, Insert );
00538         Colored += LevelIndices.size();
00539   Level++;
00540 
00541   RowMap.Comm().SumAll( &Colored, &GlobalColored, 1 );
00542   if( verbosity_ > 2 ) std::cout << "Num Globally Colored: " << GlobalColored << " from Num Global Boundary Nodes: " << BoundarySize << std::endl;
00543       }
00544     }
00545 
00546     if( verbosity_ > 1 ) std::cout << "BoundaryColoring:\n " << BoundaryColoring;
00547 
00548     Epetra_MapColoring RowColorMap( RowMap );
00549 
00550     //Add Boundary Colors
00551     for( int i = 0; i < LocalBoundarySize; ++i )
00552     {
00553       int GID = BoundaryMap.GID(i);
00554       RowColorMap(GID) = BoundaryColoring(GID);
00555     }
00556 
00557     Epetra_MapColoring Adj2ColColorMap( Adj2->ColMap() );
00558     Epetra_Import Adj2Import( Adj2->ColMap(), RowMap );
00559     Adj2ColColorMap.Import( RowColorMap, Adj2Import, Insert );
00560 
00561     if( verbosity_ > 1 ) std::cout << "RowColoringMap:\n " << RowColorMap;
00562     if( verbosity_ > 1 ) std::cout << "Adj2ColColorMap:\n " << Adj2ColColorMap;
00563 
00564     std::vector<int> rowOrder( nRows );
00565     if( reordering_ == 0 || reordering_ == 1 ) 
00566     {
00567       std::multimap<int,int> adjMap;
00568       typedef std::multimap<int,int>::value_type adjMapValueType;
00569       for( int i = 0; i < nRows; ++i )
00570         adjMap.insert( adjMapValueType( Adj2->NumMyIndices(i), i ) );
00571       std::multimap<int,int>::iterator iter = adjMap.begin();
00572       std::multimap<int,int>::iterator end = adjMap.end();
00573       if( reordering_ == 0 ) //largest first (less colors)
00574       {
00575         for( int i = 1; iter != end; ++iter, ++i )
00576           rowOrder[nRows-i] = (*iter).second;
00577       }
00578       else                  //smallest first (better balance)
00579       {
00580         for( int i = 0; iter != end; ++iter, ++i )
00581           rowOrder[i] = (*iter).second;
00582       }
00583     }
00584     else if( reordering_ == 2 ) //random
00585     {
00586       for( int i = 0; i < nRows; ++i )
00587         rowOrder[i] = i;
00588 #ifdef TFLOP
00589       random_shuffle( rowOrder.begin(), rowOrder.end() );
00590 #endif
00591     }
00592 
00593     //Constrained greedy coloring of interior
00594     set<int> InteriorColors;
00595     for( int row = 0; row < nRows; ++row )
00596     {
00597       if( !RowColorMap[ rowOrder[row] ] )
00598       {
00599         Adj2->ExtractMyRowView( rowOrder[row], NumIndices, Indices );
00600 
00601         set<int> usedColors;
00602         int color;
00603         for( int i = 0; i < NumIndices; ++i )
00604         {
00605           color = Adj2ColColorMap[ Indices[i] ];
00606           if( color > 0 ) usedColors.insert( color );
00607           color = 0;
00608           int testcolor = 1;
00609           while( !color )
00610           {
00611             if( !usedColors.count( testcolor ) ) color = testcolor;
00612             ++testcolor;
00613           }
00614         }
00615         Adj2ColColorMap[ rowOrder[row] ] = color;
00616   InteriorColors.insert( color );
00617       }
00618     }
00619     if( verbosity_ > 1 ) std::cout << MyPID << " Num Interior Colors: " << InteriorColors.size() << std::endl;
00620     if( verbosity_ > 1 ) std::cout << "RowColorMap after Greedy:\n " << RowColorMap;
00621 
00622     ColorMap = new Epetra_MapColoring( ColMap );
00623     Epetra_Import ColImport( ColMap, Adj2->ColMap() );
00624     ColorMap->Import( Adj2ColColorMap, ColImport, Insert );
00625 
00626     if( !distance1_ ) delete Adj2;
00627   }
00628 
00629   if( verbosity_ > 0 ) std::cout << MyPID << " ColorMap Color Count: " << ColorMap->NumColors() << std::endl;
00630   if( verbosity_ > 1 ) std::cout << "ColorMap!\n" << *ColorMap;
00631 
00632   newObj_ = ColorMap;
00633 
00634   return *ColorMap;
00635 }
00636 
00637 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines