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);
102 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') {
112 int MyPID = Comm.
MyPID();
115 if (verbose && MyPID==0)
118 if (verbose) cout << Comm <<endl;
123 verbose = verbose && (MyPID==0);
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);
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;
185 if ((i<nx-1) && (j>0 )) Indices[k++] = i+1 +(j-1)*nx;
193 if ((i<nx-2) && (j>0 )) Indices[k++] = i+2 +(j-1)*nx;
207 if (i<nx-1) Indices[k++] = i+1 + j *nx;
208 if (j<ny-1) Indices[k++] = i +(j+1)*nx;
216 if ((i>0 ) && (j<ny-1)) Indices[k++] = i-1 +(j+1)*nx;
223 if ((i>1 ) && (j<ny-1)) Indices[k++] = i-2 +(j+1)*nx;
232 if (verbose) cout <<
"\n\nCompleting A" << endl<< endl;
234 if (verbose) cout <<
"\n\nCompleting L0" << endl<< endl;
236 if (verbose) cout <<
"\n\nCompleting U0" << endl<< endl;
238 if (verbose) cout <<
"\n\nCompleting L1" << endl<< endl;
240 if (verbose) cout <<
"\n\nCompleting U1" << endl<< endl;
242 if (verbose) cout <<
"\n\nCompleting L2" << endl<< endl;
244 if (verbose) cout <<
"\n\nCompleting U2" << endl<< endl;
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);
323 std::vector<int> Indices1(2);
326 for (i=0; i<NumPoints1; i++)
333 else if (i == NumPoints1-1)
335 Indices1[0] = NumPoints1-1;
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;
397 assert(NumIndices==NumIndices1);
398 for (j=0; j<NumIndices1; j++) {
399 if (debug &&(Indices[j]!=Indices1[j])) {
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);
412 assert(NumIndices==NumIndices1);
413 for (j=0; j<NumIndices1; j++) {
414 if (debug &&(Indices[j]!=Indices1[j])) {
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;
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;
453 if (verbose) cout <<
"\n\nLU check OK" << endl<< endl;
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
int NumMyRows() const
Returns the number of local matrix rows.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int NumMyNonzeros() const
Returns the number of nonzero entries in the local graph.
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
virtual int MyPID() const =0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
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 NumGlobalNonzeros() const
const Epetra_BlockMap & RowMap() const
int main(int argc, char *argv[])
bool MyGID(int GID_in) const
int NumGlobalRows() const
Returns the number of global matrix rows.
const Epetra_Comm & Comm() const
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global graph.
int NumMyNonzeros() const
virtual int ConstructFilledGraph()
Does the actual construction of the graph.