MOAB: Mesh Oriented datABase  (version 5.2.1)
HypreSolver.cpp
Go to the documentation of this file.
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.org.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011 
00012 #include "HypreSolver.hpp"
00013 #include "HypreParMatrix.hpp"
00014 
00015 #ifdef MOAB_HAVE_MPI
00016 
00017 #include "hypre_parcsr.hpp"
00018 
00019 #include <fstream>
00020 #include <iostream>
00021 #include <cmath>
00022 #include <cstdlib>
00023 #include <cassert>
00024 
00025 using namespace std;
00026 
00027 namespace moab
00028 {
00029 #define moab_hypre_assert( a, b ) \
00030     {                             \
00031     }
00032 #define moab_hypre_assert_t( a, b )                         \
00033     {                                                       \
00034         if( !a )                                            \
00035         {                                                   \
00036             std::cout << "HYPRE Error: " << b << std::endl; \
00037             exit( -1 );                                     \
00038         }
00039 
00040 HypreSolver::HypreSolver( bool iterative ) : AbstractSolver( iterative )
00041 {
00042     A            = NULL;
00043     setup_called = 0;
00044     B = X = NULL;
00045 }
00046 
00047 HypreSolver::HypreSolver( HypreParMatrix* _A, bool iterative ) : AbstractSolver( iterative )
00048 {
00049     A            = _A;
00050     setup_called = 0;
00051     B = X = NULL;
00052 }
00053 
00054 HYPRE_Int HypreSolver::Solve( const HypreParVector& b, HypreParVector& x ) const
00055 {
00056     HYPRE_Int err;
00057 
00058     if( A == NULL ) { MB_SET_ERR_RET_VAL( "HypreSolver::Solve (...) : HypreParMatrix A is missing", 1 ); }
00059 
00060     if( !setup_called )
00061     {
00062         err          = SetupFcn()( *this, *A, b.x_par, x.x_par );
00063         setup_called = 1;
00064     }
00065 
00066     if( !iterative_mode ) { x = 0.0; }
00067 
00068     err = SolveFcn()( *this, *A, b.x_par, x.x_par );
00069     return err;
00070 }
00071 
00072 HypreSolver::~HypreSolver()
00073 {
00074     if( B ) { delete B; }
00075 
00076     if( X ) { delete X; }
00077 }
00078 
00079 HyprePCG::HyprePCG( HypreParMatrix& _A ) : HypreSolver( &_A, true )
00080 {
00081     iterative_mode = true;
00082     HYPRE_ParCSRPCGCreate( A->GetParallelCommunicator()->comm(), &pcg_solver );
00083 }
00084 
00085 void HyprePCG::SetTol( double tol )
00086 {
00087     HYPRE_PCGSetTol( pcg_solver, tol );
00088 }
00089 
00090 void HyprePCG::SetMaxIter( int max_iter )
00091 {
00092     HYPRE_PCGSetMaxIter( pcg_solver, max_iter );
00093 }
00094 
00095 void HyprePCG::SetLogging( int logging )
00096 {
00097     HYPRE_PCGSetLogging( pcg_solver, logging );
00098 }
00099 
00100 void HyprePCG::Verbosity( int print_lvl )
00101 {
00102     HYPRE_ParCSRPCGSetPrintLevel( pcg_solver, print_lvl );
00103 }
00104 
00105 void HyprePCG::SetPreconditioner( HypreSolver& precond )
00106 {
00107     HYPRE_ParCSRPCGSetPrecond( pcg_solver, precond.SolveFcn(), precond.SetupFcn(), precond );
00108 }
00109 
00110 void HyprePCG::SetResidualConvergenceOptions( int res_frequency, double rtol )
00111 {
00112     HYPRE_PCGSetTwoNorm( pcg_solver, 1 );
00113 
00114     if( res_frequency > 0 ) { HYPRE_PCGSetRecomputeResidualP( pcg_solver, res_frequency ); }
00115 
00116     if( rtol > 0.0 ) { HYPRE_PCGSetResidualTol( pcg_solver, rtol ); }
00117 }
00118 
00119 HYPRE_Int HyprePCG::Solve( const HypreParVector& b, HypreParVector& x ) const
00120 {
00121     HYPRE_Int err;
00122     int myid;
00123     HYPRE_Int time_index = 0;
00124     HYPRE_Int num_iterations;
00125     double final_res_norm;
00126     moab::ParallelComm* pcomm;
00127     HYPRE_Int print_level;
00128     HYPRE_PCGGetPrintLevel( pcg_solver, &print_level );
00129     pcomm = A->GetParallelCommunicator();
00130 
00131     if( !setup_called )
00132     {
00133         if( print_level > 0 )
00134         {
00135 #ifdef HYPRE_TIMING
00136             time_index = hypre_InitializeTiming( "PCG Setup" );
00137             hypre_BeginTiming( time_index );
00138 #endif
00139         }
00140 
00141         err          = HYPRE_ParCSRPCGSetup( pcg_solver, *A, (HYPRE_ParVector)b, (HYPRE_ParVector)x );
00142         setup_called = 1;
00143 
00144         if( err != 0 )
00145         {
00146             cout << "PCG Setup failed with error = " << err << endl;
00147             return err;
00148         }
00149 
00150         if( print_level > 0 )
00151         {
00152 #ifdef HYPRE_TIMING
00153             hypre_EndTiming( time_index );
00154             hypre_PrintTiming( "Setup phase times", pcomm->comm() );
00155             hypre_FinalizeTiming( time_index );
00156             hypre_ClearTiming();
00157 #endif
00158         }
00159     }
00160 
00161     if( print_level > 0 )
00162     {
00163 #ifdef HYPRE_TIMING
00164         time_index = hypre_InitializeTiming( "PCG Solve" );
00165         hypre_BeginTiming( time_index );
00166 #endif
00167     }
00168 
00169     if( !iterative_mode ) { x = 0.0; }
00170 
00171     err = HYPRE_ParCSRPCGSolve( pcg_solver, *A, (HYPRE_ParVector)b, (HYPRE_ParVector)x );
00172 
00173     if( print_level > 0 )
00174     {
00175 #ifdef HYPRE_TIMING
00176         hypre_EndTiming( time_index );
00177         hypre_PrintTiming( "Solve phase times", pcomm->comm() );
00178         hypre_FinalizeTiming( time_index );
00179         hypre_ClearTiming();
00180 #endif
00181         HYPRE_ParCSRPCGGetNumIterations( pcg_solver, &num_iterations );
00182         HYPRE_ParCSRPCGGetFinalRelativeResidualNorm( pcg_solver, &final_res_norm );
00183         MPI_Comm_rank( pcomm->comm(), &myid );
00184 
00185         if( myid == 0 )
00186         {
00187             cout << "PCG Iterations = " << num_iterations << endl
00188                  << "Final PCG Relative Residual Norm = " << final_res_norm << endl;
00189         }
00190     }
00191 
00192     return err;
00193 }
00194 
00195 HyprePCG::~HyprePCG()
00196 {
00197     HYPRE_ParCSRPCGDestroy( pcg_solver );
00198 }
00199 
00200 HypreGMRES::HypreGMRES( HypreParMatrix& _A ) : HypreSolver( &_A, true )
00201 {
00202     MPI_Comm comm;
00203     int k_dim      = 50;
00204     int max_iter   = 100;
00205     double tol     = 1e-6;
00206     iterative_mode = true;
00207     HYPRE_ParCSRMatrixGetComm( *A, &comm );
00208     HYPRE_ParCSRGMRESCreate( comm, &gmres_solver );
00209     HYPRE_ParCSRGMRESSetKDim( gmres_solver, k_dim );
00210     HYPRE_ParCSRGMRESSetMaxIter( gmres_solver, max_iter );
00211     HYPRE_ParCSRGMRESSetTol( gmres_solver, tol );
00212 }
00213 
00214 void HypreGMRES::SetTol( double atol, double rtol )
00215 {
00216     HYPRE_GMRESSetAbsoluteTol( gmres_solver, atol );
00217     HYPRE_GMRESSetTol( gmres_solver, rtol );
00218 }
00219 
00220 void HypreGMRES::SetMaxIter( int max_iter )
00221 {
00222     HYPRE_GMRESSetMaxIter( gmres_solver, max_iter );
00223 }
00224 
00225 void HypreGMRES::SetKDim( int k_dim )
00226 {
00227     HYPRE_GMRESSetKDim( gmres_solver, k_dim );
00228 }
00229 
00230 void HypreGMRES::SetLogging( int logging )
00231 {
00232     HYPRE_GMRESSetLogging( gmres_solver, logging );
00233 }
00234 
00235 void HypreGMRES::Verbosity( int print_lvl )
00236 {
00237     HYPRE_GMRESSetPrintLevel( gmres_solver, print_lvl );
00238 }
00239 
00240 void HypreGMRES::SetPreconditioner( HypreSolver& precond )
00241 {
00242     HYPRE_ParCSRGMRESSetPrecond( gmres_solver, precond.SolveFcn(), precond.SetupFcn(), precond );
00243 }
00244 
00245 HYPRE_Int HypreGMRES::Solve( const HypreParVector& b, HypreParVector& x ) const
00246 {
00247     HYPRE_Int err;
00248     int myid;
00249     HYPRE_Int time_index = 0;
00250     HYPRE_Int num_iterations;
00251     double final_res_norm;
00252     MPI_Comm comm;
00253     HYPRE_Int print_level;
00254     HYPRE_GMRESGetPrintLevel( gmres_solver, &print_level );
00255     HYPRE_ParCSRMatrixGetComm( *A, &comm );
00256 
00257     if( !setup_called )
00258     {
00259         if( print_level > 0 )
00260         {
00261 #ifdef HYPRE_TIMING
00262             time_index = hypre_InitializeTiming( "GMRES Setup" );
00263             hypre_BeginTiming( time_index );
00264 #endif
00265         }
00266 
00267         err          = HYPRE_ParCSRGMRESSetup( gmres_solver, *A, b, x );
00268         setup_called = 1;
00269 
00270         if( err != 0 )
00271         {
00272             cout << "PCG Setup failed with error = " << err << endl;
00273             return err;
00274         }
00275 
00276         if( print_level > 0 )
00277         {
00278 #ifdef HYPRE_TIMING
00279             hypre_EndTiming( time_index );
00280             hypre_PrintTiming( "Setup phase times", comm );
00281             hypre_FinalizeTiming( time_index );
00282             hypre_ClearTiming();
00283 #endif
00284         }
00285     }
00286 
00287     if( print_level > 0 )
00288     {
00289 #ifdef HYPRE_TIMING
00290         time_index = hypre_InitializeTiming( "GMRES Solve" );
00291         hypre_BeginTiming( time_index );
00292 #endif
00293     }
00294 
00295     if( !iterative_mode ) { x = 0.0; }
00296 
00297     err = HYPRE_ParCSRGMRESSolve( gmres_solver, *A, b, x );
00298 
00299     if( print_level > 0 )
00300     {
00301 #ifdef HYPRE_TIMING
00302         hypre_EndTiming( time_index );
00303         hypre_PrintTiming( "Solve phase times", comm );
00304         hypre_FinalizeTiming( time_index );
00305         hypre_ClearTiming();
00306 #endif
00307         HYPRE_ParCSRGMRESGetNumIterations( gmres_solver, &num_iterations );
00308         HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm( gmres_solver, &final_res_norm );
00309         MPI_Comm_rank( comm, &myid );
00310 
00311         if( myid == 0 )
00312         {
00313             cout << "GMRES Iterations = " << num_iterations << endl
00314                  << "Final GMRES Relative Residual Norm = " << final_res_norm << endl;
00315         }
00316     }
00317 
00318     return err;
00319 }
00320 
00321 HypreGMRES::~HypreGMRES()
00322 {
00323     HYPRE_ParCSRGMRESDestroy( gmres_solver );
00324 }
00325 
00326 HypreParaSails::HypreParaSails( HypreParMatrix& A ) : HypreSolver( &A )
00327 {
00328     MPI_Comm comm;
00329     int sai_max_levels   = 1;
00330     double sai_threshold = 0.1;
00331     double sai_filter    = 0.1;
00332     int sai_sym          = 0;
00333     double sai_loadbal   = 0.0;
00334     int sai_reuse        = 0;
00335     int sai_logging      = 1;
00336     HYPRE_ParCSRMatrixGetComm( A, &comm );
00337     HYPRE_ParaSailsCreate( comm, &sai_precond );
00338     HYPRE_ParaSailsSetParams( sai_precond, sai_threshold, sai_max_levels );
00339     HYPRE_ParaSailsSetFilter( sai_precond, sai_filter );
00340     HYPRE_ParaSailsSetSym( sai_precond, sai_sym );
00341     HYPRE_ParaSailsSetLoadbal( sai_precond, sai_loadbal );
00342     HYPRE_ParaSailsSetReuse( sai_precond, sai_reuse );
00343     HYPRE_ParaSailsSetLogging( sai_precond, sai_logging );
00344 }
00345 
00346 void HypreParaSails::SetSymmetry( int sym )
00347 {
00348     HYPRE_ParaSailsSetSym( sai_precond, sym );
00349 }
00350 
00351 HypreParaSails::~HypreParaSails()
00352 {
00353     HYPRE_ParaSailsDestroy( sai_precond );
00354 }
00355 
00356 HypreBoomerAMG::HypreBoomerAMG()
00357 {
00358     HYPRE_BoomerAMGCreate( &amg_precond );
00359     SetDefaultOptions();
00360 }
00361 
00362 HypreBoomerAMG::HypreBoomerAMG( HypreParMatrix& A ) : HypreSolver( &A )
00363 {
00364     HYPRE_BoomerAMGCreate( &amg_precond );
00365     SetDefaultOptions();
00366 }
00367 
00368 void HypreBoomerAMG::SetDefaultOptions()
00369 {
00370     // AMG coarsening options:
00371     int coarsen_type = 10;    // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
00372     int agg_levels   = 2;     // number of aggressive coarsening levels
00373     double theta     = 0.25;  // strength threshold: 0.25, 0.5, 0.8
00374     // AMG interpolation options:
00375     int interp_type = 6;  // 6 = extended+i, 0 = classical
00376     int Pmax        = 4;  // max number of elements per row in P
00377     // AMG relaxation options:
00378     int relax_type   = 10;  // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
00379     int relax_sweeps = 1;   // relaxation sweeps on each level
00380     // Additional options:
00381     int print_level = 1;   // print AMG iterations? 1 = no, 2 = yes
00382     int max_levels  = 25;  // max number of levels in AMG hierarchy
00383     HYPRE_BoomerAMGSetCoarsenType( amg_precond, coarsen_type );
00384     HYPRE_BoomerAMGSetAggNumLevels( amg_precond, agg_levels );
00385     HYPRE_BoomerAMGSetRelaxType( amg_precond, relax_type );
00386     HYPRE_BoomerAMGSetNumSweeps( amg_precond, relax_sweeps );
00387     HYPRE_BoomerAMGSetStrongThreshold( amg_precond, theta );
00388     HYPRE_BoomerAMGSetInterpType( amg_precond, interp_type );
00389     HYPRE_BoomerAMGSetPMaxElmts( amg_precond, Pmax );
00390     HYPRE_BoomerAMGSetPrintLevel( amg_precond, print_level );
00391     HYPRE_BoomerAMGSetMaxLevels( amg_precond, max_levels );
00392     // Use as a preconditioner (one V-cycle, zero tolerance)
00393     HYPRE_BoomerAMGSetMaxIter( amg_precond, 1 );
00394     HYPRE_BoomerAMGSetTol( amg_precond, 0.0 );
00395 }
00396 
00397 void HypreBoomerAMG::ResetAMGPrecond()
00398 {
00399     HYPRE_Int coarsen_type;
00400     HYPRE_Int agg_levels;
00401     HYPRE_Int relax_type;
00402     HYPRE_Int relax_sweeps;
00403     HYPRE_Real theta;
00404     HYPRE_Int interp_type;
00405     HYPRE_Int Pmax;
00406     HYPRE_Int print_level;
00407     HYPRE_Int dim;
00408     hypre_ParAMGData* amg_data = (hypre_ParAMGData*)amg_precond;
00409     // read options from amg_precond
00410     HYPRE_BoomerAMGGetCoarsenType( amg_precond, &coarsen_type );
00411     agg_levels   = hypre_ParAMGDataAggNumLevels( amg_data );
00412     relax_type   = hypre_ParAMGDataUserRelaxType( amg_data );
00413     relax_sweeps = hypre_ParAMGDataUserNumSweeps( amg_data );
00414     HYPRE_BoomerAMGGetStrongThreshold( amg_precond, &theta );
00415     hypre_BoomerAMGGetInterpType( amg_precond, &interp_type );
00416     HYPRE_BoomerAMGGetPMaxElmts( amg_precond, &Pmax );
00417     HYPRE_BoomerAMGGetPrintLevel( amg_precond, &print_level );
00418     HYPRE_BoomerAMGGetNumFunctions( amg_precond, &dim );
00419     HYPRE_BoomerAMGDestroy( amg_precond );
00420     HYPRE_BoomerAMGCreate( &amg_precond );
00421     HYPRE_BoomerAMGSetCoarsenType( amg_precond, coarsen_type );
00422     HYPRE_BoomerAMGSetAggNumLevels( amg_precond, agg_levels );
00423     HYPRE_BoomerAMGSetRelaxType( amg_precond, relax_type );
00424     HYPRE_BoomerAMGSetNumSweeps( amg_precond, relax_sweeps );
00425     HYPRE_BoomerAMGSetMaxLevels( amg_precond, 25 );
00426     HYPRE_BoomerAMGSetTol( amg_precond, 0.0 );
00427     HYPRE_BoomerAMGSetMaxIter( amg_precond, 1 );  // one V-cycle
00428     HYPRE_BoomerAMGSetStrongThreshold( amg_precond, theta );
00429     HYPRE_BoomerAMGSetInterpType( amg_precond, interp_type );
00430     HYPRE_BoomerAMGSetPMaxElmts( amg_precond, Pmax );
00431     HYPRE_BoomerAMGSetPrintLevel( amg_precond, print_level );
00432     HYPRE_BoomerAMGSetNumFunctions( amg_precond, dim );
00433 }
00434 
00435 void HypreBoomerAMG::SetOperator( const HypreParMatrix& op )
00436 {
00437     if( A ) { ResetAMGPrecond(); }
00438 
00439     // update base classes: Operator, Solver, HypreSolver
00440     A            = const_cast< HypreParMatrix* >( &op );
00441     setup_called = 0;
00442     delete X;
00443     delete B;
00444     B = X = NULL;
00445 }
00446 
00447 void HypreBoomerAMG::SetSystemsOptions( int dim )
00448 {
00449     HYPRE_BoomerAMGSetNumFunctions( amg_precond, dim );
00450     // More robust options with respect to convergence
00451     HYPRE_BoomerAMGSetAggNumLevels( amg_precond, 0 );
00452     HYPRE_BoomerAMGSetStrongThreshold( amg_precond, dim * 0.25 );
00453 }
00454 
00455 HypreBoomerAMG::~HypreBoomerAMG()
00456 {
00457     HYPRE_BoomerAMGDestroy( amg_precond );
00458 }
00459 
00460 }  // namespace moab
00461 
00462 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines