60int 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 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.
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[])