71 int main(
int argc,
char *argv[])
83 MPI_Init(&argc,&argv);
102 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') {
112 int MyPID = Comm.
MyPID();
118 if (
verbose) cout << Comm <<endl;
125 if (verbose1 && argc != 4) {
128 cout <<
"Setting nx = " << nx <<
", ny = " << ny << endl;
130 else if (!verbose1 && argc != 3) {
135 nx = atoi(argv[nextarg++]);
136 if (nx<3) {cout <<
"nx = " << nx <<
": Must be greater than 2 for meaningful graph." << endl; exit(1);}
137 ny = atoi(argv[nextarg++]);
138 if (ny<3) {cout <<
"ny = " << ny <<
": Must be greater than 2 for meaningful graph." << endl; exit(1);}
141 long long NumGlobalPoints = nx*ny;
142 long long IndexBase = 0;
145 cout <<
"\n\n*****Building 5 point matrix, Level 1 and 2 filled matrices for" << endl
146 <<
" nx = " << nx <<
", ny = " << ny << endl<< endl;
151 Epetra_Map Map(NumGlobalPoints, IndexBase, Comm);
165 std::vector<long long> Indices(4);
167 for (j=0; j<ny; j++) {
168 for (i=0; i<nx; i++) {
169 long long Row = i+j*nx;
170 if (Map.
MyGID(Row)) {
177 if (i>0) Indices[k++] = i-1 + j *nx;
178 if (j>0) Indices[k++] = i +(j-1)*nx;
182 assert(L0.InsertGlobalIndices(Row, k, &Indices[0])>=0);
185 if ((i<nx-1) && (j>0 )) Indices[k++] = i+1 +(j-1)*nx;
189 assert(L1.InsertGlobalIndices(Row, k, &Indices[0])>=0);
193 if ((i<nx-2) && (j>0 )) Indices[k++] = i+2 +(j-1)*nx;
196 assert(L2.InsertGlobalIndices(Row, k, &Indices[0])>=0);
207 if (i<nx-1) Indices[k++] = i+1 + j *nx;
208 if (j<ny-1) Indices[k++] = i +(j+1)*nx;
212 assert(U0.InsertGlobalIndices(Row, k, &Indices[0])>=0);
216 if ((i>0 ) && (j<ny-1)) Indices[k++] = i-1 +(j+1)*nx;
219 assert(U1.InsertGlobalIndices(Row, k, &Indices[0])>=0);
223 if ((i>1 ) && (j<ny-1)) Indices[k++] = i-2 +(j+1)*nx;
226 assert(U2.InsertGlobalIndices(Row, k, &Indices[0])>=0);
232 if (
verbose) cout <<
"\n\nCompleting A" << endl<< endl;
234 if (
verbose) cout <<
"\n\nCompleting L0" << endl<< endl;
235 assert(L0.FillComplete()==0);
236 if (
verbose) cout <<
"\n\nCompleting U0" << endl<< endl;
237 assert(U0.FillComplete()==0);
238 if (
verbose) cout <<
"\n\nCompleting L1" << endl<< endl;
239 assert(L1.FillComplete()==0);
240 if (
verbose) cout <<
"\n\nCompleting U1" << endl<< endl;
241 assert(U1.FillComplete()==0);
242 if (
verbose) cout <<
"\n\nCompleting L2" << endl<< endl;
243 assert(L2.FillComplete()==0);
244 if (
verbose) cout <<
"\n\nCompleting U2" << endl<< endl;
245 assert(U2.FillComplete()==0);
247 if (
verbose) cout <<
"\n\n*****Testing ILU(0) constructor on A" << endl<< endl;
250 assert(ILU0.ConstructFilledGraph()==0);
252 assert(
check(L0, U0, ILU0, NumGlobalPoints, NumMyPoints, 0,
verbose)==0);
254 if (
verbose) cout <<
"\n\n*****Testing ILU(1) constructor on A" << endl<< endl;
257 assert(ILU1.ConstructFilledGraph()==0);
259 assert(
check(L1, U1, ILU1, NumGlobalPoints, NumMyPoints, 1,
verbose)==0);
261 if (
verbose) cout <<
"\n\n*****Testing ILU(2) constructor on A" << endl<< endl;
264 assert(ILU2.ConstructFilledGraph()==0);
266 assert(
check(L2, U2, ILU2, NumGlobalPoints, NumMyPoints, 2,
verbose)==0);
268 if (
verbose) cout <<
"\n\n*****Testing copy constructor" << endl<< endl;
272 assert(
check(L2, U2, ILUC, NumGlobalPoints, NumMyPoints, 2,
verbose)==0);
274 if (
verbose) cout <<
"\n\n*****Testing copy constructor" << endl<< endl;
276 Teuchos::RefCountPtr<Ifpack_IlukGraph> OverlapGraph;
277 for (
int overlap = 1; overlap < 4; overlap++) {
278 if (
verbose) cout <<
"\n\n*********************************************" << endl;
279 if (
verbose) cout <<
"\n\nConstruct Level 1 fill with Overlap = " << overlap <<
".\n\n" << endl;
282 assert(OverlapGraph->ConstructFilledGraph()==0);
285 cout <<
"Number of Global Rows = " << OverlapGraph->NumGlobalRows64() << endl;
286 cout <<
"Number of Global Nonzeros = " << OverlapGraph->NumGlobalNonzeros64() << endl;
287 cout <<
"Number of Local Rows = " << OverlapGraph->NumMyRows() << endl;
288 cout <<
"Number of Local Nonzeros = " << OverlapGraph->NumMyNonzeros() << endl;
296 int NumElements1 = 6;
297 long long NumPoints1 = NumElements1;
302 std::vector<int> NumNz1(NumPoints1);
307 for (i=0; i<NumPoints1; i++)
308 if (i==0 || i == NumPoints1-1)
315 Epetra_Map Map1(NumPoints1, NumPoints1, 1, Comm);
323 std::vector<long long> Indices1(2);
326 for (i=0; i<NumPoints1; i++)
333 else if (i == NumPoints1-1)
335 Indices1[0] = NumPoints1-1;
344 assert(A1.InsertGlobalIndices(i+1, NumEntries1, &Indices1[0])==0);
346 assert(A1.InsertGlobalIndices(ip1, 1, &ip1)==0);
350 assert(A1.FillComplete()==0);
352 if (
verbose) cout <<
"\n\nPrint out tridiagonal matrix with IndexBase = 1.\n\n" << endl;
355 if (
verbose) cout <<
"\n\nConstruct Level 1 fill with IndexBase = 1.\n\n" << endl;
358 assert(ILU11.ConstructFilledGraph()==0);
360 if (
verbose) cout <<
"\n\nPrint out Level 1 ILU graph of tridiagonal matrix with IndexBase = 1.\n\n" << endl;
361 if (verbose1) cout << ILU11 << endl;
376 long long NumGlobalRows1,
int NumMyRows1,
int LevelFill1,
bool verbose) {
381 int NumIndices, * Indices;
382 int NumIndices1, * Indices1;
393 for (i=0; i<LU.NumMyRows(); i++) {
396 assert(L1.ExtractMyRowView(i, NumIndices1, Indices1)==0);
397 assert(NumIndices==NumIndices1);
398 for (j=0; j<NumIndices1; j++) {
399 if (debug &&(Indices[j]!=Indices1[j])) {
400 int MyPID = L.
RowMap().Comm().MyPID();
401 cout <<
"Proc " << MyPID
402 <<
" Local Row = " << i
403 <<
" L.Indices["<< j <<
"] = " << Indices[j]
404 <<
" L1.Indices["<< j <<
"] = " << Indices1[j] << endl;
406 assert(Indices[j]==Indices1[j]);
408 Nout += (NumIndices-NumIndices1);
411 assert(U1.ExtractMyRowView(i, NumIndices1, Indices1)==0);
412 assert(NumIndices==NumIndices1);
413 for (j=0; j<NumIndices1; j++) {
414 if (debug &&(Indices[j]!=Indices1[j])) {
415 int MyPID = L.
RowMap().Comm().MyPID();
416 cout <<
"Proc " << MyPID
417 <<
" Local Row = " << i
418 <<
" U.Indices["<< j <<
"] = " << Indices[j]
419 <<
" U1.Indices["<< j <<
"] = " << Indices1[j] << endl;
421 assert(Indices[j]==Indices1[j]);
423 Nout += (NumIndices-NumIndices1);
428 long long NumGlobalRows = LU.NumGlobalRows64();
429 if (
verbose) cout <<
"\n\nNumber of Global Rows = " << NumGlobalRows << endl<< endl;
431 assert(NumGlobalRows==NumGlobalRows1);
433 long long NumGlobalNonzeros = LU.NumGlobalNonzeros64();
434 if (
verbose) cout <<
"\n\nNumber of Global Nonzero entries = "
435 << NumGlobalNonzeros << endl<< endl;
439 L.
RowMap().Comm().SumAll(&Nout, &NoutG, 1);
443 int NumMyRows = LU.NumMyRows();
444 if (
verbose) cout <<
"\n\nNumber of Rows = " << NumMyRows << endl<< endl;
446 assert(NumMyRows==NumMyRows1);
448 int NumMyNonzeros = LU.NumMyNonzeros();
449 if (
verbose) cout <<
"\n\nNumber of Nonzero entries = " << NumMyNonzeros << endl<< endl;
453 if (
verbose) cout <<
"\n\nLU check OK" << endl<< endl;
int main(int argc, char *argv[])
int check(Epetra_CrsGraph &L, Epetra_CrsGraph &U, Ifpack_IlukGraph &LU, long long NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose)