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 long long 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();
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;
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->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;
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 long long 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;
std::string Ifpack_Version()
const Epetra_Comm & Comm() const
bool MyGID(int GID_in) const
virtual int MyPID() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
const Epetra_BlockMap & RowMap() const
long long NumGlobalNonzeros64() const
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
int NumMyNonzeros() const
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
long long NumGlobalNonzeros64() const
Returns the number of nonzero entries in the global graph.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
long long NumGlobalRows64() const
Returns the number of global matrix rows.
int NumMyRows() const
Returns the number of local matrix rows.
int NumMyNonzeros() const
Returns the number of nonzero entries in the local graph.
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
virtual int ConstructFilledGraph()
Does the actual construction of the graph.
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)