MOAB: Mesh Oriented datABase  (version 5.3.1)
hypre_test.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines