Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/IlukGraph/cxx_main.cpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43// FIXME: assert below are not correct if NDEBUG is defined
44#ifdef NDEBUG
45#undef NDEBUG
46#endif
47
48#include "Ifpack_ConfigDefs.h"
49#include <stdio.h>
50#include <stdlib.h>
51#include <ctype.h>
52#include <assert.h>
53#include <string.h>
54#include <math.h>
55#include "Epetra_Map.h"
56#include "Ifpack_Version.h"
57#include "Epetra_CrsGraph.h"
58#include "Ifpack_IlukGraph.h"
59#include "Teuchos_RefCountPtr.hpp"
60#ifdef EPETRA_MPI
61#include "Epetra_MpiComm.h"
62#else
63#include "Epetra_SerialComm.h"
64#endif
65
66// Prototype
67
69 int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose);
70
71 int main(int argc, char *argv[])
72{
73 using std::cout;
74 using std::endl;
75
76 int ierr = 0, i, j;
77 int nx, ny;
78
79#ifdef EPETRA_MPI
80
81 // Initialize MPI
82
83 MPI_Init(&argc,&argv);
84 //int size, rank; // Number of MPI processes, My process ID
85
86 //MPI_Comm_size(MPI_COMM_WORLD, &size);
87 //MPI_Comm_rank(MPI_COMM_WORLD, &rank);
88 Epetra_MpiComm Comm( MPI_COMM_WORLD );
89
90#else
91
92 //int size = 1; // Serial case (not using MPI)
93 //int rank = 0;
94
96#endif
97
98 bool verbose = false;
99
100 int nextarg = 1;
101 // Check if we should print results to standard out
102 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') {
103 verbose = true;
104 nextarg++;
105 }
106
107 // char tmp;
108 // if (rank==0) cout << "Press any key to continue..."<< endl;
109 // if (rank==0) cin >> tmp;
110 // Comm.Barrier();
111
112 int MyPID = Comm.MyPID();
113 int NumProc = Comm.NumProc();
114
115 if (verbose && MyPID==0)
116 cout << Ifpack_Version() << endl << endl;
117
118 if (verbose) cout << Comm <<endl;
119
120 //int sqrtNumProc = (int) ceil(sqrt((double) NumProc));
121
122 bool verbose1 = verbose;
123 verbose = verbose && (MyPID==0);
124
125 if (verbose1 && argc != 4) {
126 nx = 10;
127 ny = 12*NumProc;
128 cout << "Setting nx = " << nx << ", ny = " << ny << endl;
129 }
130 else if (!verbose1 && argc != 3) {
131 nx = 10;
132 ny = 12*NumProc;
133 }
134 else {
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);}
139 }
140
141 int NumGlobalPoints = nx*ny;
142 int IndexBase = 0;
143
144 if (verbose)
145 cout << "\n\n*****Building 5 point matrix, Level 1 and 2 filled matrices for" << endl
146 << " nx = " << nx << ", ny = " << ny << endl<< endl;
147
148
149 // Create a 5 point stencil graph, level 1 fill of it and level 2 fill of it
150
151 Epetra_Map Map(NumGlobalPoints, IndexBase, Comm);
152
153 int NumMyPoints = Map.NumMyPoints();
154
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);
162
163 // Add rows one-at-a-time
164
165 std::vector<int> Indices(4); // Work space
166
167 for (j=0; j<ny; j++) {
168 for (i=0; i<nx; i++) {
169 int Row = i+j*nx;
170 if (Map.MyGID(Row)) { // Only work on rows I own
171
172 //**** Work on lower triangle of all three matrices ****
173
174 // Define entries (i-1,j), (i,j-1)
175
176 int k = 0;
177 if (i>0) Indices[k++] = i-1 + j *nx;
178 if (j>0) Indices[k++] = i +(j-1)*nx;
179
180 // Define lower triangular terms of original matrix and L(0)
181 assert(A.InsertGlobalIndices(Row, k, &Indices[0])>=0);
182 assert(L0.InsertGlobalIndices(Row, k, &Indices[0])>=0);
183
184 // Define entry (i+1,j-1)
185 if ((i<nx-1) && (j>0 )) Indices[k++] = i+1 +(j-1)*nx;
186
187
188 // Define lower triangle of level(1) fill matrix
189 assert(L1.InsertGlobalIndices(Row, k, &Indices[0])>=0);
190
191 // Define entry (i+2, j-1)
192
193 if ((i<nx-2) && (j>0 )) Indices[k++] = i+2 +(j-1)*nx;
194
195 // Define lower triangle of level(2) fill matrix
196 assert(L2.InsertGlobalIndices(Row, k, &Indices[0])>=0);
197
198 // Define main diagonal of original matrix
199 assert(A.InsertGlobalIndices(Row, 1, &Row)>=0);
200
201 k = 0; // Reset index counter
202
203 //**** Work on upper triangle of all three matrices ****
204
205 // Define entries (i+1,j), ( i,j+1)
206
207 if (i<nx-1) Indices[k++] = i+1 + j *nx;
208 if (j<ny-1) Indices[k++] = i +(j+1)*nx;
209
210 // Define upper triangular terms of original matrix and L(0)
211 assert(A.InsertGlobalIndices(Row, k, &Indices[0])>=0);
212 assert(U0.InsertGlobalIndices(Row, k, &Indices[0])>=0);
213
214 // Define entry (i-1,j+1)
215
216 if ((i>0 ) && (j<ny-1)) Indices[k++] = i-1 +(j+1)*nx;
217
218 // Define upper triangle of level(1) fill matrix
219 assert(U1.InsertGlobalIndices(Row, k, &Indices[0])>=0);
220
221 // Define entry (i-2, j+1)
222
223 if ((i>1 ) && (j<ny-1)) Indices[k++] = i-2 +(j+1)*nx;
224
225 // Define upper triangle of level(2) fill matrix
226 assert(U2.InsertGlobalIndices(Row, k, &Indices[0])>=0);
227 }
228 }
229 }
230
231 // Finish up
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);
246
247 if (verbose) cout << "\n\n*****Testing ILU(0) constructor on A" << endl<< endl;
248
249 Ifpack_IlukGraph ILU0(A, 0, 0);
250 assert(ILU0.ConstructFilledGraph()==0);
251
252 assert(check(L0, U0, ILU0, NumGlobalPoints, NumMyPoints, 0, verbose)==0);
253
254 if (verbose) cout << "\n\n*****Testing ILU(1) constructor on A" << endl<< endl;
255
256 Ifpack_IlukGraph ILU1(A, 1, 0);
257 assert(ILU1.ConstructFilledGraph()==0);
258
259 assert(check(L1, U1, ILU1, NumGlobalPoints, NumMyPoints, 1, verbose)==0);
260
261 if (verbose) cout << "\n\n*****Testing ILU(2) constructor on A" << endl<< endl;
262
263 Ifpack_IlukGraph ILU2(A, 2, 0);
264 assert(ILU2.ConstructFilledGraph()==0);
265
266 assert(check(L2, U2, ILU2, NumGlobalPoints, NumMyPoints, 2, verbose)==0);
267
268 if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
269
270 Ifpack_IlukGraph ILUC(ILU2);
271
272 assert(check(L2, U2, ILUC, NumGlobalPoints, NumMyPoints, 2, verbose)==0);
273
274 if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
275
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;
280
281 OverlapGraph = Teuchos::rcp( new Ifpack_IlukGraph(A, 1, overlap) );
282 assert(OverlapGraph->ConstructFilledGraph()==0);
283
284 if (verbose) {
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;
289 }
290 }
291
292 if (verbose1) {
293 // Test ostream << operator (if verbose1)
294 // Construct a Map that puts 6 equations on each PE
295
296 int NumElements1 = 6;
297 int NumPoints1 = NumElements1;
298
299 // Create an integer vector NumNz that is used to build the Petra Matrix.
300 // NumNz[i] is the Number of terms for the ith global equation on this processor
301
302 std::vector<int> NumNz1(NumPoints1);
303
304 // We are building a tridiagonal matrix where each row has (-1 2 -1)
305 // So we need 2 off-diagonal terms (except for the first and last equation)
306
307 for (i=0; i<NumPoints1; i++)
308 if (i==0 || i == NumPoints1-1)
309 NumNz1[i] = 2;
310 else
311 NumNz1[i] = 3;
312
313 // Create a Epetra_Matrix
314
315 Epetra_Map Map1(NumPoints1, NumPoints1, 1, Comm);
316 Epetra_CrsGraph A1(Copy, Map1, &NumNz1[0]);
317
318 // Add rows one-at-a-time
319 // Need some vectors to help
320 // Off diagonal Values will always be -1
321
322
323 std::vector<int> Indices1(2);
324 int NumEntries1;
325
326 for (i=0; i<NumPoints1; i++)
327 {
328 if (i==0)
329 {
330 Indices1[0] = 2;
331 NumEntries1 = 1;
332 }
333 else if (i == NumPoints1-1)
334 {
335 Indices1[0] = NumPoints1-1;
336 NumEntries1 = 1;
337 }
338 else
339 {
340 Indices1[0] = i;
341 Indices1[1] = i+2;
342 NumEntries1 = 2;
343 }
344 assert(A1.InsertGlobalIndices(i+1, NumEntries1, &Indices1[0])==0);
345 int ip1 = i+1;
346 assert(A1.InsertGlobalIndices(ip1, 1, &ip1)==0); // Put in the diagonal entry
347 }
348
349 // Finish up
350 assert(A1.FillComplete()==0);
351
352 if (verbose) cout << "\n\nPrint out tridiagonal matrix with IndexBase = 1.\n\n" << endl;
353 cout << A1 << endl;
354
355 if (verbose) cout << "\n\nConstruct Level 1 fill with IndexBase = 1.\n\n" << endl;
356
357 Ifpack_IlukGraph ILU11(A1, 1, 0);
358 assert(ILU11.ConstructFilledGraph()==0);
359
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;
362
363 }
364
365#ifdef EPETRA_MPI
366 MPI_Finalize() ;
367#endif
368
369
370/* end main
371*/
372return ierr ;
373}
374
376 int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose) {
377 using std::cout;
378 using std::endl;
379
380 int i, j;
381 int NumIndices, * Indices;
382 int NumIndices1, * Indices1;
383
384 bool debug = true;
385
386 Epetra_CrsGraph& L1 = LU.L_Graph();
387 Epetra_CrsGraph& U1 = LU.U_Graph();
388
389 // Test entries and count nonzeros
390
391 int Nout = 0;
392
393 for (i=0; i<LU.NumMyRows(); i++) {
394
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;
405 }
406 assert(Indices[j]==Indices1[j]);
407 }
408 Nout += (NumIndices-NumIndices1);
409
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;
420 }
421 assert(Indices[j]==Indices1[j]);
422 }
423 Nout += (NumIndices-NumIndices1);
424 }
425
426 // Test query functions
427
428 int NumGlobalRows = LU.NumGlobalRows();
429 if (verbose) cout << "\n\nNumber of Global Rows = " << NumGlobalRows << endl<< endl;
430
431 assert(NumGlobalRows==NumGlobalRows1);
432
433 int NumGlobalNonzeros = LU.NumGlobalNonzeros();
434 if (verbose) cout << "\n\nNumber of Global Nonzero entries = "
435 << NumGlobalNonzeros << endl<< endl;
436
437 int NoutG = 0;
438
439 L.RowMap().Comm().SumAll(&Nout, &NoutG, 1);
440
441 assert(NumGlobalNonzeros==L.NumGlobalNonzeros()+U.NumGlobalNonzeros()-NoutG);
442
443 int NumMyRows = LU.NumMyRows();
444 if (verbose) cout << "\n\nNumber of Rows = " << NumMyRows << endl<< endl;
445
446 assert(NumMyRows==NumMyRows1);
447
448 int NumMyNonzeros = LU.NumMyNonzeros();
449 if (verbose) cout << "\n\nNumber of Nonzero entries = " << NumMyNonzeros << endl<< endl;
450
451 assert(NumMyNonzeros==L.NumMyNonzeros()+U.NumMyNonzeros()-Nout);
452
453 if (verbose) cout << "\n\nLU check OK" << endl<< endl;
454
455 return(0);
456}
Copy
std::string Ifpack_Version()
int NumMyPoints() const
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
int NumGlobalNonzeros() const
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
int NumMyNonzeros() const
int NumProc() const
int MyPID() const
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global graph.
int NumGlobalRows() const
Returns the number of global matrix rows.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
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, int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose)