MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 // Parts of this test have been taken from MFEM and HPCG benchmark example 00002 00003 #include "moab/Core.hpp" 00004 #include "HypreSolver.hpp" 00005 00006 #include <iostream> 00007 using namespace std; 00008 00009 #undef DEBUG 00010 00011 moab::ErrorCode GenerateTestMatrixAndVectors( int nx, int ny, int nz, moab::HypreParMatrix& A, 00012 moab::HypreParVector& sol, moab::HypreParVector& rhs, 00013 moab::HypreParVector& exactsol ); 00014 00015 int main( int argc, char* argv[] ) 00016 { 00017 // double norm, d; 00018 int ierr = 0; 00019 double times[5]; 00020 int nx = 5, ny = 5, nz = 5; 00021 #ifdef MOAB_HAVE_MPI 00022 MPI_Init( &argc, &argv ); 00023 int size, rank; // Number of MPI processes, My process ID 00024 MPI_Comm_size( MPI_COMM_WORLD, &size ); 00025 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00026 #else 00027 int size = 1; // Serial case (not using MPI) 00028 int rank = 0; 00029 #endif 00030 MPI_Comm comm = MPI_COMM_WORLD; 00031 00032 moab::Core mbCore; 00033 moab::Interface& mb = mbCore; 00034 00035 // rval = mb.load_file(example.c_str(), &euler_set, opts.c_str()); 00036 00037 moab::ParallelComm* pcomm = new moab::ParallelComm( &mb, comm ); 00038 00039 moab::HypreParMatrix A( pcomm ); 00040 moab::HypreParVector x( pcomm ), b( pcomm ), xexact( pcomm ); 00041 #ifdef DEBUG 00042 00043 if( rank == 0 ) 00044 { 00045 int junk = 0; 00046 cout << "Press enter to continue" << endl; 00047 cin >> junk; 00048 } 00049 00050 MPI_Barrier( MPI_COMM_WORLD ); 00051 #endif 00052 00053 if( argc > 4 ) 00054 { 00055 if( rank == 0 ) 00056 cerr << "Usage:" << endl 00057 << "\t" << argv[0] << " nx ny nz" << endl 00058 << " where nx, ny and nz are the local sub-block dimensions, or" << endl; 00059 00060 exit( 1 ); 00061 } 00062 else 00063 { 00064 if( rank == 0 ) 00065 cout << "Using command:" 00066 << "\t" << argv[0] << " nx ny nz" << endl; 00067 } 00068 00069 times[0] = MPI_Wtime(); // Initialize it (if needed) 00070 00071 if( argc == 4 ) 00072 { 00073 nx = atoi( argv[1] ); 00074 ny = atoi( argv[2] ); 00075 nz = atoi( argv[3] ); 00076 } 00077 00078 GenerateTestMatrixAndVectors( nx, ny, nz, A, x, b, xexact ); 00079 00080 times[1] = MPI_Wtime() - times[0]; // operator creation time 00081 const bool useCG = true; // TODO make into command line option 00082 const bool useAMG = false; // TODO make into command line option 00083 int niters = 0; 00084 double normr = 0; 00085 int max_iter = 150; 00086 double tolerance = 1e-15; // Set tolerance to zero to make all runs do max_iter iterations 00087 times[2] = MPI_Wtime(); // reset 00088 moab::HypreSolver *lsSolver, *psSolver; 00089 00090 if( useCG ) 00091 { 00092 lsSolver = new moab::HyprePCG( A ); 00093 moab::HyprePCG* cgSolver = dynamic_cast< moab::HyprePCG* >( lsSolver ); 00094 cgSolver->SetTol( tolerance ); 00095 cgSolver->SetMaxIter( max_iter ); 00096 cgSolver->SetLogging( 1 ); 00097 cgSolver->Verbosity( 1 ); 00098 cgSolver->SetResidualConvergenceOptions( 3, 0.0 ); 00099 } 00100 else 00101 { 00102 lsSolver = new moab::HypreGMRES( A ); 00103 moab::HypreGMRES* gmresSolver = dynamic_cast< moab::HypreGMRES* >( lsSolver ); 00104 gmresSolver->SetTol( tolerance, normr ); 00105 gmresSolver->SetMaxIter( max_iter ); 00106 gmresSolver->SetLogging( 1 ); 00107 gmresSolver->Verbosity( 5 ); 00108 gmresSolver->SetKDim( 30 ); 00109 } 00110 00111 if( useAMG ) 00112 { 00113 psSolver = new moab::HypreBoomerAMG( A ); 00114 dynamic_cast< moab::HypreBoomerAMG* >( psSolver )->SetSystemsOptions( 3 ); 00115 } 00116 else 00117 { 00118 psSolver = new moab::HypreParaSails( A ); 00119 dynamic_cast< moab::HypreParaSails* >( psSolver )->SetSymmetry( 1 ); 00120 } 00121 00122 if( psSolver ) lsSolver->SetPreconditioner( *psSolver ); 00123 00124 times[2] = MPI_Wtime() - times[2]; // solver setup time 00125 /* Now let us solve the linear system */ 00126 times[3] = MPI_Wtime(); // reset 00127 ierr = lsSolver->Solve( b, x ); 00128 00129 if( ierr ) cerr << "Error in call to CG/GMRES HYPRE Solver: " << ierr << ".\n" << endl; 00130 00131 times[3] = MPI_Wtime() - times[3]; // solver time 00132 lsSolver->GetNumIterations( niters ); 00133 lsSolver->GetFinalResidualNorm( normr ); 00134 times[4] = MPI_Wtime() - times[0]; 00135 // x.Print ( "solution.out" ); 00136 #ifdef MOAB_HAVE_MPI 00137 double t4 = times[3]; 00138 double t4min = 0.0; 00139 double t4max = 0.0; 00140 double t4avg = 0.0; 00141 MPI_Allreduce( &t4, &t4min, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD ); 00142 MPI_Allreduce( &t4, &t4max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD ); 00143 MPI_Allreduce( &t4, &t4avg, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); 00144 t4avg = t4avg / ( (double)size ); 00145 #endif 00146 00147 if( rank == 0 ) 00148 { // Only PE 0 needs to compute and report timing results 00149 // double fniters = niters; 00150 // double fnrow = A.NNZ(); 00151 // double fnnz = A.M(); 00152 // double fnops_ddot = fniters * 4 * fnrow; 00153 // double fnops_waxpby = fniters * 6 * fnrow; 00154 // double fnops_sparsemv = fniters * 2 * fnnz; 00155 // double fnops = fnops_ddot + fnops_waxpby + fnops_sparsemv; 00156 std::cout << "Testing Laplace problem solver with HYPRE interface\n"; 00157 { 00158 std::cout << "\tParallelism\n"; 00159 #ifdef MOAB_HAVE_MPI 00160 std::cout << "Number of MPI ranks: \t" << size << std::endl; 00161 #else 00162 std::cout << "MPI not enabled\n"; 00163 #endif 00164 } 00165 std::cout << "Dimensions (nx*ny*nz): \t(" << nx << "*" << ny << "*" << nz << ")\n"; 00166 std::cout << "Number of iterations: \t" << niters << std::endl; 00167 std::cout << "Final residual: \t" << normr << std::endl; 00168 std::cout << "************ Performance Summary (times in sec) *************" << std::endl; 00169 { 00170 std::cout << "\nTime Summary" << std::endl; 00171 std::cout << "Total \t" << times[4] << std::endl; 00172 std::cout << "Operator Creation \t" << times[1] << std::endl; 00173 std::cout << "Solver Setup \t" << times[2] << std::endl; 00174 std::cout << "HYPRE Solver \t" << times[3] << std::endl; 00175 std::cout << "Avg. Solver \t" << t4avg << std::endl; 00176 } 00177 } 00178 00179 // Compute difference between known exact solution and computed solution 00180 // All processors are needed here. 00181 // double residual = 0; 00182 // if ((ierr = ComputeLinfError(x, xexact, &residual))) 00183 // cerr << "Error in call to ComputeLinfError: " << ierr << ".\n" << endl; 00184 // if (rank==0) 00185 // cout << "Difference between computed and exact = " 00186 // << residual << ".\n" << endl; 00187 delete lsSolver; 00188 // Finish up 00189 #ifdef MOAB_HAVE_MPI 00190 MPI_Finalize(); 00191 #endif 00192 return 0; 00193 } 00194 00195 moab::ErrorCode GenerateTestMatrixAndVectors( int nx, int ny, int nz, moab::HypreParMatrix& A, 00196 moab::HypreParVector& sol, moab::HypreParVector& rhs, 00197 moab::HypreParVector& exactsol ) 00198 { 00199 #ifdef DEBUG 00200 int debug = 1; 00201 #else 00202 int debug = 0; 00203 #endif 00204 #ifdef MOAB_HAVE_MPI 00205 int size, rank; // Number of MPI processes, My process ID 00206 MPI_Comm_size( MPI_COMM_WORLD, &size ); 00207 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00208 #else 00209 int size = 1; // Serial case (not using MPI) 00210 int rank = 0; 00211 #endif 00212 // Set this bool to true if you want a 7-pt stencil instead of a 27 pt stencil 00213 bool use_7pt_stencil = false; 00214 const int NNZPERROW = 27; 00215 00216 int local_nrow = nx * ny * nz; // This is the size of our subblock 00217 assert( local_nrow > 0 ); // Must have something to work with 00218 int local_nnz = NNZPERROW * local_nrow; // Approximately 27 nonzeros per row (except for boundary nodes) 00219 int total_nrow = 0; 00220 MPI_Allreduce( &local_nrow, &total_nrow, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); 00221 00222 // int total_nrow = local_nrow * size; // Total number of grid points in mesh 00223 long long total_nnz = 00224 NNZPERROW * (long long)total_nrow; // Approximately 27 nonzeros per row (except for boundary nodes) 00225 int start_row = local_nrow * rank; // Each processor gets a section of a chimney stack domain 00226 int stop_row = start_row + local_nrow - 1; 00227 00228 // Allocate arrays that are of length local_nrow 00229 // Eigen::SparseMatrix<double, Eigen::RowMajor> diagMatrix(local_nrow, local_nrow); 00230 // diagMatrix.reserve(local_nnz); 00231 HYPRE_Int col[2] = { start_row, stop_row }; 00232 A.resize( total_nrow, col ); 00233 00234 sol.resize( total_nrow, start_row, stop_row ); 00235 rhs.resize( total_nrow, start_row, stop_row ); 00236 exactsol.resize( total_nrow, start_row, stop_row ); 00237 long long nnzglobal = 0; 00238 00239 for( int iz = 0; iz < nz; iz++ ) 00240 { 00241 for( int iy = 0; iy < ny; iy++ ) 00242 { 00243 for( int ix = 0; ix < nx; ix++ ) 00244 { 00245 const int curlocalrow = iz * nx * ny + iy * nx + ix; 00246 const int currow = start_row + curlocalrow; 00247 int nnzrow = 0; 00248 std::vector< double > colvals; 00249 std::vector< int > indices; 00250 00251 for( int sz = -1; sz <= 1; sz++ ) 00252 { 00253 for( int sy = -1; sy <= 1; sy++ ) 00254 { 00255 for( int sx = -1; sx <= 1; sx++ ) 00256 { 00257 const int curlocalcol = sz * nx * ny + sy * nx + sx; 00258 const int curcol = currow + curlocalcol; 00259 00260 // Since we have a stack of nx by ny by nz domains , stacking in the z 00261 // direction, we check to see if sx and sy are reaching outside of the 00262 // domain, while the check for the curcol being valid is sufficient to 00263 // check the z values 00264 if( ( ix + sx >= 0 ) && ( ix + sx < nx ) && ( iy + sy >= 0 ) && ( iy + sy < ny ) && 00265 ( curcol >= 0 && curcol < total_nrow ) ) 00266 { 00267 if( !use_7pt_stencil || ( sz * sz + sy * sy + sx * sx <= 1 ) ) 00268 { // This logic will skip over point that are not part of a 7-pt 00269 // stencil 00270 if( curcol == currow ) { colvals.push_back( 27.0 ); } 00271 else 00272 { 00273 colvals.push_back( -1.0 ); 00274 } 00275 00276 indices.push_back( curcol ); 00277 nnzrow++; 00278 } 00279 } 00280 } // end sx loop 00281 } // end sy loop 00282 } // end sz loop 00283 00284 int ncols = indices.size(); 00285 A.SetValues( 1, &ncols, &currow, indices.data(), colvals.data() ); 00286 // diagMatrix.set_row_values ( curlocalrow, indices.size(), indices.data(), 00287 // colvals.data() ); 00288 nnzglobal += nnzrow; 00289 sol.SetValue( currow, 0.0 ); 00290 rhs.SetValue( currow, 27.0 - ( (double)( nnzrow - 1 ) ) ); 00291 exactsol.SetValue( currow, 1.0 ); 00292 } // end ix loop 00293 } // end iy loop 00294 } // end iz loop 00295 00296 if( debug ) 00297 cout << "Global size of the matrix: " << total_nrow << " and NNZ (estimate) = " << total_nnz 00298 << " NNZ (actual) = " << nnzglobal << endl; 00299 00300 if( debug ) cout << "Process " << rank << " of " << size << " has " << local_nrow << endl; 00301 00302 if( debug ) cout << " rows. Global rows " << start_row << " through " << stop_row << endl; 00303 00304 if( debug ) cout << "Process " << rank << " of " << size << " has " << local_nnz << " nonzeros." << endl; 00305 00306 A.FinalizeAssembly(); 00307 sol.FinalizeAssembly(); 00308 rhs.FinalizeAssembly(); 00309 exactsol.FinalizeAssembly(); 00310 // A.Print("matrix.out"); 00311 return moab::MB_SUCCESS; 00312 }