Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/IlukGraph_LL/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 long long 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 long long NumGlobalPoints = nx*ny;
142 long long 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<long long> Indices(4); // Work space
166
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)) { // 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->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;
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 long long 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<long long> 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 long long 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 long long 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 long long NumGlobalRows = LU.NumGlobalRows64();
429 if (verbose) cout << "\n\nNumber of Global Rows = " << NumGlobalRows << endl<< endl;
430
431 assert(NumGlobalRows==NumGlobalRows1);
432
433 long long NumGlobalNonzeros = LU.NumGlobalNonzeros64();
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.NumGlobalNonzeros64()+U.NumGlobalNonzeros64()-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}
457
Copy
std::string Ifpack_Version()
int NumMyPoints() const
bool MyGID(int GID_in) const
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
int NumProc() const
int MyPID() const
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
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)