Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
example/MapColoring/cxx_main.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Epetra: Linear Algebra Services Package
5 // Copyright 2011 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #include <assert.h>
43 #include <stdio.h>
44 #include <set>
45 #include <algorithm>
46 #include <vector>
47 #include "Epetra_Map.h"
48 #include "Epetra_BlockMap.h"
49 #include "Epetra_MapColoring.h"
50 #include "Epetra_CrsMatrix.h"
51 #ifdef EPETRA_MPI
52 #include "mpi.h"
53 #include "Epetra_MpiComm.h"
54 #else
55 #include "Epetra_SerialComm.h"
56 #endif
57 
58 using namespace std;
59 
60 int main(int argc, char * argv[]) {
61 
62 // Set up the Epetra communicator
63 #ifdef EPETRA_MPI
64  MPI_Init(&argc, &argv);
65  int size; // Number of MPI processes
66  int rank; // My process ID
67  MPI_Comm_size(MPI_COMM_WORLD, &size);
68  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
69  Epetra_MpiComm Comm(MPI_COMM_WORLD);
70 #else
71  int size = 1; // Serial case (not using MPI)
72  int rank = 0; // Processor ID
73  Epetra_SerialComm Comm;
74 #endif
75  // cout << Comm << endl;
76 
77  int MyPID = Comm.MyPID();
78  int NumProc = Comm.NumProc();
79  bool verbose = (MyPID == 0);
80  cout << "MyPID = " << MyPID << endl
81  << "NumProc = " << NumProc << endl;
82 
83  // Get the problem size from the command line argument
84  if (argc < 2 || argc > 3) {
85  if (verbose)
86  cout << "Usage: " << argv[0] << " nx [ny]" << endl;
87  exit(1);
88  } // end if
89  int nx = atoi(argv[1]); // Get the dimensions for a 1D or 2D
90  int ny = 1; // central difference problem
91  if (argc == 3) ny = atoi(argv[2]);
92  int NumGlobalElements = nx * ny;
93  if (NumGlobalElements < NumProc) {
94  if (verbose)
95  cout << "numGlobalElements = " << NumGlobalElements
96  << " cannot be < number of processors = " << NumProc << endl;
97  exit(1);
98  } // end if
99 
100  // Epetra distribution map
101  int IndexBase = 0; // Zero-based indices
102  Epetra_Map Map(NumGlobalElements, IndexBase, Comm);
103  // if (verbose) cout << Map << endl;
104 
105  // Extract the global indices of the elements local to this processor
106  int NumMyElements = Map.NumMyElements();
107  int * MyGlobalElements = new int[NumMyElements];
108  Map.MyGlobalElements(MyGlobalElements);
109  for (int p = 0; p < NumProc; p++)
110  if (p == MyPID) {
111  cout << endl << "Processor " << MyPID << ": Global elements = ";
112  for (int i = 0; i < NumMyElements; i++)
113  cout << MyGlobalElements[i] << " ";
114  cout << endl;
115  Comm.Barrier();
116  } // end if
117 
118  // Create the number of non-zeros for a tridiagonal (1D problem) or banded
119  // (2D problem) matrix
120  int * NumNz = new int[NumMyElements];
121  int global_i;
122  int global_j;
123  for (int i = 0; i < NumMyElements; i++) {
124  NumNz[i] = 5;
125  global_j = MyGlobalElements[i] / nx;
126  global_i = MyGlobalElements[i] - global_j * nx;
127  if (global_i == 0) NumNz[i] -= 1; // By having separate statements,
128  if (global_i == nx-1) NumNz[i] -= 1; // this works for 2D as well as 1D
129  if (global_j == 0) NumNz[i] -= 1; // systems (i.e. nx x 1 or 1 x ny)
130  if (global_j == ny-1) NumNz[i] -= 1; // or even a 1 x 1 system
131  }
132  if (verbose) {
133  cout << endl << "NumNz: ";
134  for (int i = 0; i < NumMyElements; i++) cout << NumNz[i] << " ";
135  cout << endl;
136  } // end if
137 
138  // Create the Epetra Compressed Row Sparse matrix
139  // Note: the actual values for the matrix entries are meaningless for
140  // this exercise, but I assign them anyway.
141  Epetra_CrsMatrix A(Copy, Map, NumNz);
142 
143  double * Values = new double[4];
144  for (int i = 0; i < 4; i++) Values[i] = -1.0;
145  int * Indices = new int[4];
146  double diag = 2.0;
147  if (ny > 1) diag = 4.0;
148  int NumEntries;
149 
150  for (int i = 0; i < NumMyElements; i++) {
151  global_j = MyGlobalElements[i] / nx;
152  global_i = MyGlobalElements[i] - global_j * nx;
153  NumEntries = 0;
154  if (global_j > 0 && ny > 1)
155  Indices[NumEntries++] = global_i + (global_j-1)*nx;
156  if (global_i > 0)
157  Indices[NumEntries++] = global_i-1 + global_j *nx;
158  if (global_i < nx-1)
159  Indices[NumEntries++] = global_i+1 + global_j *nx;
160  if (global_j < ny-1 && ny > 1)
161  Indices[NumEntries++] = global_i + (global_j+1)*nx;
162 
163  // Put in the off-diagonal terms
164  assert(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values,
165  Indices) == 0);
166  // Put in the diagonal entry
167  assert(A.InsertGlobalValues(MyGlobalElements[i], 1, &diag,
168  MyGlobalElements+i) == 0);
169  } // end i loop
170 
171  // Finish up matrix construction
172  delete [] Values;
173  delete [] Indices;
174  assert(A.FillComplete() == 0);
175  // cout << endl << A << endl;
176 
177  // Create the local distance-1 adjancency graph
178  // This is essentially a transpose of the Epetra_CrsGraph, where off-
179  // processor couplings are ignored and global indexes are converted to
180  // local. We use the C++ standard libraries vector and set, since we
181  // don't know how many nonzeroes we will end up with for each column.
182 
183  vector< set<int> > adj1(NumMyElements);
184  for (int lr = 0; lr < adj1.size(); lr++) {
185  int lrid; // Local row ID
186  double * Values = new double[NumNz[lr]];
187  int * Indices = new int[NumNz[lr]];
188  assert(A.ExtractMyRowCopy(lr, NumNz[lr], NumNz[lr], Values, Indices) == 0);
189  for (int i = 0; i < NumNz[lr]; i++) {
190  lrid = A.LRID(Indices[i]);
191  if (lrid >= 0) adj1[lrid].insert(lr);
192  } // end i loop
193  delete [] Values;
194  delete [] Indices;
195  } // end lr loop
196 
197  if (verbose) {
198  cout << endl;
199  for (int lr = 0; lr < NumMyElements; lr++) {
200  cout << "adj1[" << lr << "] = { ";
201  for (set<int>::const_iterator p = adj1[lr].begin(); p != adj1[lr].end();
202  p++) cout << *p << " ";
203  cout << "}" << endl;
204  } // end lr loop
205  } // end if
206 
207  // Create the local distance-2 adjancency graph
208  // This is built from the distance-1 adjancency graph. We use the C++
209  // standard libraries vector and set, since we don't know how many
210  // nonzeroes we will end up with for each column.
211 
212  vector< set<int> > adj2(NumMyElements);
213  for (int lc = 0; lc < NumMyElements; lc++) {
214  for (set<int>::const_iterator p = adj1[lc].begin(); p != adj1[lc].end();
215  p++) {
216  int lrid; // Local row ID
217  double * Values = new double[NumNz[*p]];
218  int * Indices = new int[NumNz[*p]];
219  assert(A.ExtractMyRowCopy(*p, NumNz[*p], NumNz[*p], Values, Indices)
220  == 0);
221  for (int i = 0; i < NumNz[*p]; i++) {
222  lrid = A.LRID(Indices[i]);
223  if (lrid >= 0) adj2[lc].insert(lrid);
224  } // end i loop
225  delete [] Values;
226  delete [] Indices;
227  } // end p loop
228  } // end lc loop
229 
230  cout << endl;
231  for (int lc = 0; lc < NumMyElements; lc++) {
232  cout << "adj2[" << lc << "] = { ";
233  for (set<int>::const_iterator p = adj2[lc].begin(); p != adj2[lc].end();
234  p++) cout << *p << " ";
235  cout << "}" << endl;
236  } // end lc loop
237 
238  // Now that we have the local distance-2 adjacency graph, we can compute a
239  // color map using a greedy algorithm. The first step is to compute Delta,
240  // the maximum size (degree) of adj1.
241  size_t Delta = 0;
242  for (int i = 0; i < NumMyElements; i++)
243  Delta = max(Delta, adj1[i].size());
244  cout << endl << "Delta = " << Delta << endl << endl;
245 
246  // Now create a color map and initialize all values to 0, which
247  // indicates that none of the columns have yet been colored.
248  int * color_map = new int[NumMyElements];
249  for (int i = 0; i < NumMyElements; i++) color_map[i] = 0;
250 
251  // Apply the distance-2 greedy coloring algorithm
252  for (int column = 0; column < NumMyElements; column++) {
253  set<int> allowedColors; // Create the set of allowed colors
254  for (int i = 1; i < Delta*Delta+1; i++) allowedColors.insert(i);
255  for (set<int>::const_iterator p = adj2[column].begin();
256  p != adj2[column].end(); p++) if (color_map[*p] > 0)
257  allowedColors.erase(color_map[*p]);
258  color_map[column] = *(allowedColors.begin());
259  cout << "color_map[" << column << "] = " << color_map[column] << endl;
260  } // end col loop
261 
262  // New section to Epetra_MapColoring
263  Epetra_MapColoring C1(Map, color_map);
264 
265  cout << C1;
266 
267  // Clean up
268 
269  delete [] MyGlobalElements;
270  delete [] NumNz;
271  delete [] color_map;
272 
273  cout << endl << argv[0] << " done." << endl;
274 
275 } // end main
276 
int LRID(int GRID_in) const
Returns the local row index for given global row index, returns -1 if no local row for this global ro...
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_MapColoring: A class for coloring Epetra_Map and Epetra_BlockMap objects.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given global row of the matrix.
int MyPID() const
Return my process ID.
Epetra_MpiComm: The Epetra MPI Communication Class.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
int NumMyElements() const
Number of elements on the calling processor.
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
Epetra_SerialComm: The Epetra Serial Communication Class.
void Barrier() const
Epetra_SerialComm Barrier function.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int main(int argc, char *argv[])
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.