55 #include "Epetra_Map.h" 57 #include "Epetra_CrsGraph.h" 59 #include "Teuchos_RefCountPtr.hpp" 61 #include "Epetra_MpiComm.h" 63 #include "Epetra_SerialComm.h" 69 int NumGlobalRows1,
int NumMyRows1,
int LevelFill1,
bool verbose);
71 int main(
int argc,
char *argv[])
83 MPI_Init(&argc,&argv);
88 Epetra_MpiComm Comm( MPI_COMM_WORLD );
95 Epetra_SerialComm Comm;
102 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') {
112 int MyPID = Comm.MyPID();
113 int NumProc = Comm.NumProc();
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 int NumGlobalPoints = nx*ny;
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);
153 int NumMyPoints = Map.NumMyPoints();
155 Epetra_CrsGraph
A(
Copy, Map, 5);
156 Epetra_CrsGraph L0(
Copy, Map, Map, 2);
157 Epetra_CrsGraph U0(
Copy, Map, Map, 2);
158 Epetra_CrsGraph L1(
Copy, Map, Map, 3);
159 Epetra_CrsGraph U1(
Copy, Map, Map, 3);
160 Epetra_CrsGraph L2(
Copy, Map, Map, 4);
161 Epetra_CrsGraph U2(
Copy, Map, Map, 4);
165 std::vector<int> Indices(4);
167 for (j=0; j<ny; j++) {
168 for (i=0; i<nx; i++) {
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;
181 assert(
A.InsertGlobalIndices(Row, k, &Indices[0])>=0);
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);
199 assert(
A.InsertGlobalIndices(Row, 1, &Row)>=0);
207 if (i<nx-1) Indices[k++] = i+1 + j *nx;
208 if (j<ny-1) Indices[k++] = i +(j+1)*nx;
211 assert(
A.InsertGlobalIndices(Row, k, &Indices[0])>=0);
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;
233 assert(
A.FillComplete()==0);
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;
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;
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;
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->NumGlobalRows() << endl;
286 cout <<
"Number of Global Nonzeros = " << OverlapGraph->NumGlobalNonzeros() << endl;
287 cout <<
"Number of Local Rows = " << OverlapGraph->NumMyRows() << endl;
288 cout <<
"Number of Local Nonzeros = " << OverlapGraph->NumMyNonzeros() << endl;
296 int NumElements1 = 6;
297 int 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);
316 Epetra_CrsGraph A1(
Copy, Map1, &NumNz1[0]);
323 std::vector<int> 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;
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 int NumGlobalRows1,
int NumMyRows1,
int LevelFill1,
bool verbose) {
381 int NumIndices, * Indices;
382 int NumIndices1, * Indices1;
386 Epetra_CrsGraph& L1 = LU.
L_Graph();
387 Epetra_CrsGraph& U1 = LU.
U_Graph();
395 assert(L.ExtractMyRowView(i, NumIndices, Indices)==0);
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);
410 assert(U.ExtractMyRowView(i, NumIndices, Indices)==0);
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);
429 if (
verbose) cout <<
"\n\nNumber of Global Rows = " << NumGlobalRows << endl<< endl;
431 assert(NumGlobalRows==NumGlobalRows1);
434 if (
verbose) cout <<
"\n\nNumber of Global Nonzero entries = " 435 << NumGlobalNonzeros << endl<< endl;
439 L.RowMap().Comm().SumAll(&Nout, &NoutG, 1);
441 assert(NumGlobalNonzeros==L.NumGlobalNonzeros()+U.NumGlobalNonzeros()-NoutG);
444 if (
verbose) cout <<
"\n\nNumber of Rows = " << NumMyRows << endl<< endl;
446 assert(NumMyRows==NumMyRows1);
449 if (
verbose) cout <<
"\n\nNumber of Nonzero entries = " << NumMyNonzeros << endl<< endl;
451 assert(NumMyNonzeros==L.NumMyNonzeros()+U.NumMyNonzeros()-Nout);
453 if (
verbose) cout <<
"\n\nLU check OK" << endl<< endl;
int NumGlobalRows() const
Returns the number of global matrix rows.
int NumMyRows() const
Returns the number of local matrix rows.
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global graph.
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int main(int argc, char *argv[])
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int check(Epetra_CrsGraph &L, Epetra_CrsGraph &U, Ifpack_IlukGraph &LU, int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose)
std::string Ifpack_Version()
int NumMyNonzeros() const
Returns the number of nonzero entries in the local graph.
virtual int ConstructFilledGraph()
Does the actual construction of the graph.