EpetraExt Package Browser (Single Doxygen Collection) Development
test/MapColoring/cxx_main.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 // CrsMatrix_MapColoring Test routine
00043 #include <Epetra_ConfigDefs.h>
00044 #include "EpetraExt_Version.h"
00045 
00046 #ifdef EPETRA_MPI
00047 #include "Epetra_MpiComm.h"
00048 #include <mpi.h>
00049 #endif
00050 
00051 #include "Epetra_SerialComm.h"
00052 #include "Epetra_Time.h"
00053 #include "Epetra_Map.h"
00054 #include "Epetra_CrsGraph.h"
00055 #include "Epetra_IntVector.h"
00056 #include "Epetra_MapColoring.h"
00057 #include "EpetraExt_MapColoring.h"
00058 #include "EpetraExt_MapColoringIndex.h"
00059 #include "../epetra_test_err.h"
00060 
00061 #include <vector>
00062 
00063 // prototype
00064 void printColoring (const Epetra_MapColoring & ColorMap, Epetra_CrsGraph *A, bool verbose);
00065 
00066 int main(int argc, char *argv[]) {
00067 
00068 
00069 #ifdef EPETRA_MPI
00070 
00071   // Initialize MPI
00072 
00073   MPI_Init( &argc, &argv );
00074   //int size, rank; // Number of MPI processes, My process ID
00075 
00076   //MPI_Comm_size(MPI_COMM_WORLD, &size);
00077   //MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00078 
00079 #else
00080 
00081   //int size = 1; // Serial case (not using MPI)
00082   //int rank = 0;
00083 
00084 #endif
00085 
00086   bool verbose = false;
00087 
00088   int nx = 5;
00089   int ny = 5;
00090 
00091   if( argc > 1 )
00092   {
00093     if( argc > 4 )
00094     {
00095       cout << "Usage: " << argv[0] << " [-v [nx [ny]]]" << endl;
00096       exit(1);
00097     }
00098 
00099     int loc = 1;
00100     // Check if we should print results to standard out
00101     if(argv[loc][0]=='-' && argv[loc][1]=='v')
00102     { verbose = true; ++loc; }
00103 
00104     if (loc < argc) nx = atoi( argv[loc++] );
00105     if( loc < argc) ny = atoi( argv[loc] );
00106   }
00107 
00108 #ifdef EPETRA_MPI
00109   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00110 #else
00111   Epetra_SerialComm Comm;
00112 #endif
00113 
00114   int MyPID = Comm.MyPID();
00115   int NumProc = Comm.NumProc();
00116 
00117   bool verbose1 = false;
00118   if(verbose) verbose1 = (MyPID==0);
00119 
00120   if(verbose1)
00121     cout << EpetraExt::EpetraExt_Version() << endl << endl;
00122 
00123   Comm.Barrier();
00124 
00125   if(verbose) cout << Comm << endl << flush;
00126   Comm.Barrier();
00127 
00128   int NumGlobalElements = nx * ny;
00129   if( NumGlobalElements < NumProc )
00130   {
00131     cout << "NumGlobalElements = " << NumGlobalElements <<
00132             " cannot be < number of processors = " << NumProc;
00133     exit(1);
00134   } 
00135   
00136   int IndexBase = 0;
00137   Epetra_Map Map( NumGlobalElements, IndexBase, Comm );
00138 
00139   // Extract the global indices of the elements local to this processor
00140   int NumMyElements = Map.NumMyElements();
00141   std::vector<int> MyGlobalElements( NumMyElements );
00142   Map.MyGlobalElements( &MyGlobalElements[0] );
00143   if( verbose ) cout << Map;
00144 
00145   // Create the number of non-zeros for a tridiagonal (1D problem) or banded
00146   // (2D problem) matrix
00147   std::vector<int> NumNz( NumMyElements, 5 );
00148   int global_i;
00149   int global_j;
00150   for (int i = 0; i < NumMyElements; ++i)
00151   {
00152     global_j = MyGlobalElements[i] / nx;
00153     global_i = MyGlobalElements[i] - global_j * nx;
00154     if (global_i == 0)    NumNz[i] -= 1;  // By having separate statements,
00155     if (global_i == nx-1) NumNz[i] -= 1;  // this works for 2D as well as 1D
00156     if (global_j == 0)    NumNz[i] -= 1;  // systems (i.e. nx x 1 or 1 x ny)
00157     if (global_j == ny-1) NumNz[i] -= 1;  // or even a 1 x 1 system
00158   }
00159   if(verbose)
00160   { 
00161     cout << endl << "NumNz: ";
00162     for (int i = 0; i < NumMyElements; i++) cout << NumNz[i] << " ";
00163     cout << endl;
00164   } // end if
00165   
00166   // Create the Epetra Compressed Row Sparse Graph
00167   Epetra_CrsGraph A( Copy, Map, &NumNz[0] );
00168   
00169   std::vector<int> Indices(5);
00170   int NumEntries;
00171   
00172   for (int i = 0; i < NumMyElements; ++i )
00173   {
00174     global_j = MyGlobalElements[i] / nx;
00175     global_i = MyGlobalElements[i] - global_j * nx;
00176     NumEntries = 0;
00177     // (i,j-1) entry
00178     if (global_j > 0 && ny > 1)
00179       Indices[NumEntries++] = global_i   + (global_j-1)*nx;
00180     // (i-1,j) entry
00181     if (global_i > 0)
00182       Indices[NumEntries++] = global_i-1 +  global_j   *nx;
00183     // (i,j) entry
00184     Indices[NumEntries++] = MyGlobalElements[i];
00185     // (i+1,j) entry
00186     if (global_i < nx-1)
00187       Indices[NumEntries++] = global_i+1 +  global_j   *nx;
00188     // (i,j+1) entry
00189     if (global_j < ny-1 && ny > 1)
00190       Indices[NumEntries++] = global_i   + (global_j+1)*nx;
00191 
00192     // Insert the global indices
00193     A.InsertGlobalIndices( MyGlobalElements[i], NumEntries, &Indices[0] );
00194   } // end i loop
00195 
00196   // Finish up graph construction
00197   A.FillComplete();
00198 
00199   EpetraExt::CrsGraph_MapColoring
00200     Greedy0MapColoringTransform( EpetraExt::CrsGraph_MapColoring::GREEDY,
00201                      0, false, verbose );
00202   Epetra_MapColoring & Greedy0ColorMap = Greedy0MapColoringTransform( A );
00203   printColoring(Greedy0ColorMap, &A,verbose);
00204 
00205   EpetraExt::CrsGraph_MapColoring
00206     Greedy1MapColoringTransform( EpetraExt::CrsGraph_MapColoring::GREEDY,
00207                      1, false, verbose );
00208   Epetra_MapColoring & Greedy1ColorMap = Greedy1MapColoringTransform( A );
00209   printColoring(Greedy1ColorMap, &A,verbose);
00210 
00211   EpetraExt::CrsGraph_MapColoring
00212     Greedy2MapColoringTransform( EpetraExt::CrsGraph_MapColoring::GREEDY,
00213                      2, false, verbose );
00214   Epetra_MapColoring & Greedy2ColorMap = Greedy2MapColoringTransform( A );
00215   printColoring(Greedy2ColorMap, &A,verbose);
00216 
00217   EpetraExt::CrsGraph_MapColoring
00218     Lubi0MapColoringTransform( EpetraExt::CrsGraph_MapColoring::LUBY,
00219                    0, false, verbose );
00220   Epetra_MapColoring & Lubi0ColorMap = Lubi0MapColoringTransform( A );
00221   printColoring(Lubi0ColorMap, &A,verbose);
00222 
00223   EpetraExt::CrsGraph_MapColoring
00224     Lubi1MapColoringTransform( EpetraExt::CrsGraph_MapColoring::LUBY,
00225                    1, false, verbose );
00226   Epetra_MapColoring & Lubi1ColorMap = Lubi1MapColoringTransform( A );
00227   printColoring(Lubi1ColorMap, &A,verbose);
00228 
00229   EpetraExt::CrsGraph_MapColoring
00230     Lubi2MapColoringTransform( EpetraExt::CrsGraph_MapColoring::LUBY,
00231                    2, false, verbose );
00232   Epetra_MapColoring & Lubi2ColorMap = Lubi2MapColoringTransform( A );
00233   printColoring(Lubi2ColorMap, &A,verbose);
00234 
00235 #ifdef EPETRA_MPI
00236   if( verbose ) cout << "Parallel Map Coloring 1!\n";
00237   EpetraExt::CrsGraph_MapColoring
00238     Parallel1MapColoringTransform( EpetraExt::CrsGraph_MapColoring::PSEUDO_PARALLEL,
00239                        0, false, verbose );
00240   Epetra_MapColoring & Parallel1ColorMap = Parallel1MapColoringTransform( A );
00241   printColoring(Parallel1ColorMap, &A,verbose);
00242 
00243   if( verbose ) cout << "Parallel Map Coloring 2!\n";
00244   EpetraExt::CrsGraph_MapColoring
00245     Parallel2MapColoringTransform( EpetraExt::CrsGraph_MapColoring::JONES_PLASSMAN,
00246                        0, false, verbose );
00247   Epetra_MapColoring & Parallel2ColorMap = Parallel2MapColoringTransform( A );
00248   printColoring(Parallel2ColorMap, &A,verbose);
00249 #endif
00250 
00251 
00252 #ifdef EPETRA_MPI
00253   MPI_Finalize();
00254 #endif
00255 
00256   return 0;
00257 }
00258 
00259 void printColoring (const Epetra_MapColoring & ColorMap, Epetra_CrsGraph * A, bool verbose) {
00260 
00261   int NumColors = ColorMap.NumColors();
00262   int * ListOfColors = ColorMap.ListOfColors();
00263 
00264   EpetraExt::CrsGraph_MapColoringIndex MapColoringIndexTransform( ColorMap );
00265   vector<Epetra_IntVector> & ColIndices = MapColoringIndexTransform( *A );
00266 
00267   if( verbose )
00268   {
00269 
00270     cout << endl;
00271     cout << "***************************************\n";
00272     cout << "Column Indexing by Color:\n";
00273     cout << "***************************************\n";
00274     cout << endl;
00275     for( int i = 0; i < NumColors; ++i )
00276     {
00277       cout << "COLOR: " << ListOfColors[i] << endl;
00278       cout << ColIndices[i];
00279     }
00280     cout << endl;
00281   }
00282 
00283   return;
00284 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines