60 int main(
int argc,
char * argv[]) {
64 MPI_Init(&argc, &argv);
67 MPI_Comm_size(MPI_COMM_WORLD, &size);
68 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
77 int MyPID = Comm.
MyPID();
79 bool verbose = (MyPID == 0);
80 cout <<
"MyPID = " << MyPID << endl
81 <<
"NumProc = " << NumProc << endl;
84 if (argc < 2 || argc > 3) {
86 cout <<
"Usage: " << argv[0] <<
" nx [ny]" << endl;
89 int nx = atoi(argv[1]);
91 if (argc == 3) ny = atoi(argv[2]);
92 int NumGlobalElements = nx * ny;
93 if (NumGlobalElements < NumProc) {
95 cout <<
"numGlobalElements = " << NumGlobalElements
96 <<
" cannot be < number of processors = " << NumProc << endl;
102 Epetra_Map Map(NumGlobalElements, IndexBase, Comm);
107 int * MyGlobalElements =
new int[NumMyElements];
109 for (
int p = 0; p < NumProc; p++)
111 cout << endl <<
"Processor " << MyPID <<
": Global elements = ";
112 for (
int i = 0; i < NumMyElements; i++)
113 cout << MyGlobalElements[i] <<
" ";
120 int * NumNz =
new int[NumMyElements];
123 for (
int i = 0; i < NumMyElements; i++) {
125 global_j = MyGlobalElements[i] / nx;
126 global_i = MyGlobalElements[i] - global_j * nx;
127 if (global_i == 0) NumNz[i] -= 1;
128 if (global_i == nx-1) NumNz[i] -= 1;
129 if (global_j == 0) NumNz[i] -= 1;
130 if (global_j == ny-1) NumNz[i] -= 1;
133 cout << endl <<
"NumNz: ";
134 for (
int i = 0; i < NumMyElements; i++) cout << NumNz[i] <<
" ";
143 double * Values =
new double[4];
144 for (
int i = 0; i < 4; i++) Values[i] = -1.0;
145 int * Indices =
new int[4];
147 if (ny > 1) diag = 4.0;
150 for (
int i = 0; i < NumMyElements; i++) {
151 global_j = MyGlobalElements[i] / nx;
152 global_i = MyGlobalElements[i] - global_j * nx;
154 if (global_j > 0 && ny > 1)
155 Indices[NumEntries++] = global_i + (global_j-1)*nx;
157 Indices[NumEntries++] = global_i-1 + global_j *nx;
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;
168 MyGlobalElements+i) == 0);
183 vector< set<int> > adj1(NumMyElements);
184 for (
int lr = 0; lr < adj1.size(); lr++) {
186 double * Values =
new double[NumNz[lr]];
187 int * Indices =
new int[NumNz[lr]];
189 for (
int i = 0; i < NumNz[lr]; i++) {
190 lrid = A.
LRID(Indices[i]);
191 if (lrid >= 0) adj1[lrid].insert(lr);
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 <<
" ";
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();
217 double * Values =
new double[NumNz[*p]];
218 int * Indices =
new int[NumNz[*p]];
221 for (
int i = 0; i < NumNz[*p]; i++) {
222 lrid = A.
LRID(Indices[i]);
223 if (lrid >= 0) adj2[lc].insert(lrid);
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 <<
" ";
242 for (
int i = 0; i < NumMyElements; i++)
243 Delta = max(Delta, adj1[i].size());
244 cout << endl <<
"Delta = " << Delta << endl << endl;
248 int * color_map =
new int[NumMyElements];
249 for (
int i = 0; i < NumMyElements; i++) color_map[i] = 0;
252 for (
int column = 0; column < NumMyElements; column++) {
253 set<int> allowedColors;
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;
269 delete [] MyGlobalElements;
273 cout << endl << argv[0] <<
" done." << endl;
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.
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.