Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
58using namespace std;
59
60int 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
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 MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int NumMyElements() const
Number of elements on the calling processor.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
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...
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.
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.
Epetra_MapColoring: A class for coloring Epetra_Map and Epetra_BlockMap objects.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
Epetra_MpiComm: The Epetra MPI Communication Class.
void Barrier() const
Epetra_MpiComm Barrier function.
int NumProc() const
Returns total number of processes.
int MyPID() const
Return my process ID.
Epetra_SerialComm: The Epetra Serial Communication Class.
int main(int argc, char *argv[])