1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323 | // Parts of this test have been taken from MFEM and HPCG benchmark example
#include "moab/Core.hpp"
#include "HypreSolver.hpp"
#include <iostream>
using namespace std;
#undef DEBUG
moab::ErrorCode GenerateTestMatrixAndVectors( int nx,
int ny,
int nz,
moab::HypreParMatrix& A,
moab::HypreParVector& sol,
moab::HypreParVector& rhs,
moab::HypreParVector& exactsol );
int main( int argc, char* argv[] )
{
// double norm, d;
int ierr = 0;
double times[5];
int nx = 5, ny = 5, nz = 5;
#ifdef MOAB_HAVE_MPI
MPI_Init( &argc, &argv );
int size, rank; // Number of MPI processes, My process ID
MPI_Comm_size( MPI_COMM_WORLD, &size );
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
#else
int size = 1; // Serial case (not using MPI)<--- Variable 'size' is assigned a value that is never used.
int rank = 0;<--- 'rank' is assigned value '0' here.<--- 'rank' is assigned value '0' here.<--- 'rank' is assigned value '0' here.
#endif
MPI_Comm comm = MPI_COMM_WORLD;
moab::Core mbCore;
moab::Interface& mb = mbCore;
// rval = mb.load_file(example.c_str(), &euler_set, opts.c_str());
moab::ParallelComm* pcomm = new moab::ParallelComm( &mb, comm );
moab::HypreParMatrix A( pcomm );
moab::HypreParVector x( pcomm ), b( pcomm ), xexact( pcomm );
#ifdef DEBUG
if( rank == 0 )
{
int junk = 0;
cout << "Press enter to continue" << endl;
cin >> junk;
}
MPI_Barrier( MPI_COMM_WORLD );
#endif
if( argc > 4 )
{
if( rank == 0 )<--- The comparison 'rank == 0' is always true. [+]Finding the same expression on both sides of an operator is suspicious and might indicate a cut and paste or logic error. Please examine this code carefully to determine if it is correct.
cerr << "Usage:" << endl
<< "\t" << argv[0] << " nx ny nz" << endl
<< " where nx, ny and nz are the local sub-block dimensions, or" << endl;
exit( 1 );
}
else
{
if( rank == 0 )<--- The comparison 'rank == 0' is always true. [+]Finding the same expression on both sides of an operator is suspicious and might indicate a cut and paste or logic error. Please examine this code carefully to determine if it is correct.
cout << "Using command:"
<< "\t" << argv[0] << " nx ny nz" << endl;
}
times[0] = MPI_Wtime(); // Initialize it (if needed)
if( argc == 4 )
{
nx = atoi( argv[1] );
ny = atoi( argv[2] );
nz = atoi( argv[3] );
}
GenerateTestMatrixAndVectors( nx, ny, nz, A, x, b, xexact );
times[1] = MPI_Wtime() - times[0]; // operator creation time
const bool useCG = true; // TODO make into command line option
const bool useAMG = false; // TODO make into command line option
int niters = 0;
double normr = 0;
int max_iter = 150;
double tolerance = 1e-15; // Set tolerance to zero to make all runs do max_iter iterations
times[2] = MPI_Wtime(); // reset
moab::HypreSolver *lsSolver, *psSolver;
if( useCG )
{
lsSolver = new moab::HyprePCG( A );
moab::HyprePCG* cgSolver = dynamic_cast< moab::HyprePCG* >( lsSolver );
cgSolver->SetTol( tolerance );
cgSolver->SetMaxIter( max_iter );
cgSolver->SetLogging( 1 );
cgSolver->Verbosity( 1 );
cgSolver->SetResidualConvergenceOptions( 3, 0.0 );
}
else
{
lsSolver = new moab::HypreGMRES( A );
moab::HypreGMRES* gmresSolver = dynamic_cast< moab::HypreGMRES* >( lsSolver );
gmresSolver->SetTol( tolerance, normr );
gmresSolver->SetMaxIter( max_iter );
gmresSolver->SetLogging( 1 );
gmresSolver->Verbosity( 5 );
gmresSolver->SetKDim( 30 );
}
if( useAMG )
{
psSolver = new moab::HypreBoomerAMG( A );
dynamic_cast< moab::HypreBoomerAMG* >( psSolver )->SetSystemsOptions( 3 );
}
else
{
psSolver = new moab::HypreParaSails( A );
dynamic_cast< moab::HypreParaSails* >( psSolver )->SetSymmetry( 1 );
}
if( psSolver ) lsSolver->SetPreconditioner( *psSolver );
times[2] = MPI_Wtime() - times[2]; // solver setup time
/* Now let us solve the linear system */
times[3] = MPI_Wtime(); // reset
ierr = lsSolver->Solve( b, x );
if( ierr ) cerr << "Error in call to CG/GMRES HYPRE Solver: " << ierr << ".\n" << endl;
times[3] = MPI_Wtime() - times[3]; // solver time
lsSolver->GetNumIterations( niters );
lsSolver->GetFinalResidualNorm( normr );
times[4] = MPI_Wtime() - times[0];
// x.Print ( "solution.out" );
#ifdef MOAB_HAVE_MPI
double t4 = times[3];
double t4min = 0.0;
double t4max = 0.0;
double t4avg = 0.0;
MPI_Allreduce( &t4, &t4min, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
MPI_Allreduce( &t4, &t4max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
MPI_Allreduce( &t4, &t4avg, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
t4avg = t4avg / ( (double)size );
#endif
if( rank == 0 )<--- The comparison 'rank == 0' is always true. [+]Finding the same expression on both sides of an operator is suspicious and might indicate a cut and paste or logic error. Please examine this code carefully to determine if it is correct.
{ // Only PE 0 needs to compute and report timing results
// double fniters = niters;
// double fnrow = A.NNZ();
// double fnnz = A.M();
// double fnops_ddot = fniters * 4 * fnrow;
// double fnops_waxpby = fniters * 6 * fnrow;
// double fnops_sparsemv = fniters * 2 * fnnz;
// double fnops = fnops_ddot + fnops_waxpby + fnops_sparsemv;
std::cout << "Testing Laplace problem solver with HYPRE interface\n";
{
std::cout << "\tParallelism\n";
#ifdef MOAB_HAVE_MPI
std::cout << "Number of MPI ranks: \t" << size << std::endl;
#else
std::cout << "MPI not enabled\n";
#endif
}
std::cout << "Dimensions (nx*ny*nz): \t(" << nx << "*" << ny << "*" << nz << ")\n";
std::cout << "Number of iterations: \t" << niters << std::endl;
std::cout << "Final residual: \t" << normr << std::endl;
std::cout << "************ Performance Summary (times in sec) *************" << std::endl;
{
std::cout << "\nTime Summary" << std::endl;
std::cout << "Total \t" << times[4] << std::endl;
std::cout << "Operator Creation \t" << times[1] << std::endl;
std::cout << "Solver Setup \t" << times[2] << std::endl;
std::cout << "HYPRE Solver \t" << times[3] << std::endl;
std::cout << "Avg. Solver \t" << t4avg << std::endl;
}
}
// Compute difference between known exact solution and computed solution
// All processors are needed here.
// double residual = 0;
// if ((ierr = ComputeLinfError(x, xexact, &residual)))
// cerr << "Error in call to ComputeLinfError: " << ierr << ".\n" << endl;
// if (rank==0)
// cout << "Difference between computed and exact = "
// << residual << ".\n" << endl;
delete lsSolver;
// Finish up
#ifdef MOAB_HAVE_MPI
MPI_Finalize();
#endif
return 0;
}
moab::ErrorCode GenerateTestMatrixAndVectors( int nx,
int ny,
int nz,
moab::HypreParMatrix& A,
moab::HypreParVector& sol,
moab::HypreParVector& rhs,
moab::HypreParVector& exactsol )
{
#ifdef DEBUG
int debug = 1;
#else
int debug = 0;<--- Assignment 'debug=0', assigned value is 0
#endif
#ifdef MOAB_HAVE_MPI
int size, rank; // Number of MPI processes, My process ID
MPI_Comm_size( MPI_COMM_WORLD, &size );
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
#else
int size = 1; // Serial case (not using MPI)
int rank = 0;
#endif
// Set this bool to true if you want a 7-pt stencil instead of a 27 pt stencil
bool use_7pt_stencil = false;<--- Assignment 'use_7pt_stencil=false', assigned value is 0
const int NNZPERROW = 27;
int local_nrow = nx * ny * nz; // This is the size of our subblock
assert( local_nrow > 0 ); // Must have something to work with
int local_nnz = NNZPERROW * local_nrow; // Approximately 27 nonzeros per row (except for boundary nodes)
int total_nrow = 0;
MPI_Allreduce( &local_nrow, &total_nrow, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
// int total_nrow = local_nrow * size; // Total number of grid points in mesh
long long total_nnz =
NNZPERROW * (long long)total_nrow; // Approximately 27 nonzeros per row (except for boundary nodes)
int start_row = local_nrow * rank; // Each processor gets a section of a chimney stack domain
int stop_row = start_row + local_nrow - 1;
// Allocate arrays that are of length local_nrow
// Eigen::SparseMatrix<double, Eigen::RowMajor> diagMatrix(local_nrow, local_nrow);
// diagMatrix.reserve(local_nnz);
HYPRE_Int col[2] = { start_row, stop_row };
A.resize( total_nrow, col );
sol.resize( total_nrow, start_row, stop_row );
rhs.resize( total_nrow, start_row, stop_row );
exactsol.resize( total_nrow, start_row, stop_row );
long long nnzglobal = 0;
for( int iz = 0; iz < nz; iz++ )
{
for( int iy = 0; iy < ny; iy++ )
{
for( int ix = 0; ix < nx; ix++ )
{
const int curlocalrow = iz * nx * ny + iy * nx + ix;
const int currow = start_row + curlocalrow;
int nnzrow = 0;
std::vector< double > colvals;
std::vector< int > indices;
for( int sz = -1; sz <= 1; sz++ )
{
for( int sy = -1; sy <= 1; sy++ )
{
for( int sx = -1; sx <= 1; sx++ )
{
const int curlocalcol = sz * nx * ny + sy * nx + sx;
const int curcol = currow + curlocalcol;
// Since we have a stack of nx by ny by nz domains , stacking in the z
// direction, we check to see if sx and sy are reaching outside of the
// domain, while the check for the curcol being valid is sufficient to
// check the z values
if( ( ix + sx >= 0 ) && ( ix + sx < nx ) && ( iy + sy >= 0 ) && ( iy + sy < ny ) &&
( curcol >= 0 && curcol < total_nrow ) )
{
if( !use_7pt_stencil || ( sz * sz + sy * sy + sx * sx <= 1 ) )<--- Condition '!use_7pt_stencil' is always true
{ // This logic will skip over point that are not part of a 7-pt
// stencil
if( curcol == currow )
{
colvals.push_back( 27.0 );
}
else
{
colvals.push_back( -1.0 );
}
indices.push_back( curcol );
nnzrow++;
}
}
} // end sx loop
} // end sy loop
} // end sz loop
int ncols = indices.size();
A.SetValues( 1, &ncols, &currow, indices.data(), colvals.data() );
// diagMatrix.set_row_values ( curlocalrow, indices.size(), indices.data(),
// colvals.data() );
nnzglobal += nnzrow;
sol.SetValue( currow, 0.0 );
rhs.SetValue( currow, 27.0 - ( (double)( nnzrow - 1 ) ) );
exactsol.SetValue( currow, 1.0 );
} // end ix loop
} // end iy loop
} // end iz loop
if( debug )<--- Condition 'debug' is always false
cout << "Global size of the matrix: " << total_nrow << " and NNZ (estimate) = " << total_nnz
<< " NNZ (actual) = " << nnzglobal << endl;
if( debug ) cout << "Process " << rank << " of " << size << " has " << local_nrow << endl;<--- First condition
if( debug ) cout << " rows. Global rows " << start_row << " through " << stop_row << endl;<--- Second condition<--- First condition
if( debug ) cout << "Process " << rank << " of " << size << " has " << local_nnz << " nonzeros." << endl;<--- Second condition
A.FinalizeAssembly();
sol.FinalizeAssembly();
rhs.FinalizeAssembly();
exactsol.FinalizeAssembly();
// A.Print("matrix.out");
return moab::MB_SUCCESS;
}
|