hypre_test.cpp

Go to the documentation of this file.

// 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)
 int rank = 0;
#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 )
 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 )
 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 )
 { // 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;
#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;
 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 ) )
 { // 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 )
 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;
 
 if( debug ) cout << " rows. Global rows " << start_row << " through " << stop_row << endl;
 
 if( debug ) cout << "Process " << rank << " of " << size << " has " << local_nnz << " nonzeros." << endl;
 
 A.FinalizeAssembly();
 sol.FinalizeAssembly();
 rhs.FinalizeAssembly();
 exactsol.FinalizeAssembly();
 // A.Print("matrix.out");
 return moab::MB_SUCCESS;
}