Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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 )
00059     {
00060         MB_SET_ERR_RET_VAL( "HypreSolver::Solve (...) : HypreParMatrix A is missing", 1 );
00061     }
00062 
00063     if( !setup_called )
00064     {
00065         err          = SetupFcn()( *this, *A, b.x_par, x.x_par );
00066         setup_called = 1;
00067     }
00068 
00069     if( !iterative_mode )
00070     {
00071         x = 0.0;
00072     }
00073 
00074     err = SolveFcn()( *this, *A, b.x_par, x.x_par );
00075     return err;
00076 }
00077 
00078 HypreSolver::~HypreSolver()
00079 {
00080     if( B )
00081     {
00082         delete B;
00083     }
00084 
00085     if( X )
00086     {
00087         delete X;
00088     }
00089 }
00090 
00091 HyprePCG::HyprePCG( HypreParMatrix& _A ) : HypreSolver( &_A, true )
00092 {
00093     iterative_mode = true;
00094     HYPRE_ParCSRPCGCreate( A->GetParallelCommunicator()->comm(), &pcg_solver );
00095 }
00096 
00097 void HyprePCG::SetTol( double tol )
00098 {
00099     HYPRE_PCGSetTol( pcg_solver, tol );
00100 }
00101 
00102 void HyprePCG::SetMaxIter( int max_iter )
00103 {
00104     HYPRE_PCGSetMaxIter( pcg_solver, max_iter );
00105 }
00106 
00107 void HyprePCG::SetLogging( int logging )
00108 {
00109     HYPRE_PCGSetLogging( pcg_solver, logging );
00110 }
00111 
00112 void HyprePCG::Verbosity( int print_lvl )
00113 {
00114     HYPRE_ParCSRPCGSetPrintLevel( pcg_solver, print_lvl );
00115 }
00116 
00117 void HyprePCG::SetPreconditioner( HypreSolver& precond )
00118 {
00119     HYPRE_ParCSRPCGSetPrecond( pcg_solver, precond.SolveFcn(), precond.SetupFcn(), precond );
00120 }
00121 
00122 void HyprePCG::SetResidualConvergenceOptions( int res_frequency, double rtol )
00123 {
00124     HYPRE_PCGSetTwoNorm( pcg_solver, 1 );
00125 
00126     if( res_frequency > 0 )
00127     {
00128         HYPRE_PCGSetRecomputeResidualP( pcg_solver, res_frequency );
00129     }
00130 
00131     if( rtol > 0.0 )
00132     {
00133         HYPRE_PCGSetResidualTol( pcg_solver, rtol );
00134     }
00135 }
00136 
00137 HYPRE_Int HyprePCG::Solve( const HypreParVector& b, HypreParVector& x ) const
00138 {
00139     HYPRE_Int err;
00140     int myid;
00141     HYPRE_Int time_index = 0;
00142     HYPRE_Int num_iterations;
00143     double final_res_norm;
00144     moab::ParallelComm* pcomm;
00145     HYPRE_Int print_level;
00146     HYPRE_PCGGetPrintLevel( pcg_solver, &print_level );
00147     pcomm = A->GetParallelCommunicator();
00148 
00149     if( !setup_called )
00150     {
00151         if( print_level > 0 )
00152         {
00153 #ifdef HYPRE_TIMING
00154             time_index = hypre_InitializeTiming( "PCG Setup" );
00155             hypre_BeginTiming( time_index );
00156 #endif
00157         }
00158 
00159         err          = HYPRE_ParCSRPCGSetup( pcg_solver, *A, (HYPRE_ParVector)b, (HYPRE_ParVector)x );
00160         setup_called = 1;
00161 
00162         if( err != 0 )
00163         {
00164             cout << "PCG Setup failed with error = " << err << endl;
00165             return err;
00166         }
00167 
00168         if( print_level > 0 )
00169         {
00170 #ifdef HYPRE_TIMING
00171             hypre_EndTiming( time_index );
00172             hypre_PrintTiming( "Setup phase times", pcomm->comm() );
00173             hypre_FinalizeTiming( time_index );
00174             hypre_ClearTiming();
00175 #endif
00176         }
00177     }
00178 
00179     if( print_level > 0 )
00180     {
00181 #ifdef HYPRE_TIMING
00182         time_index = hypre_InitializeTiming( "PCG Solve" );
00183         hypre_BeginTiming( time_index );
00184 #endif
00185     }
00186 
00187     if( !iterative_mode )
00188     {
00189         x = 0.0;
00190     }
00191 
00192     err = HYPRE_ParCSRPCGSolve( pcg_solver, *A, (HYPRE_ParVector)b, (HYPRE_ParVector)x );
00193 
00194     if( print_level > 0 )
00195     {
00196 #ifdef HYPRE_TIMING
00197         hypre_EndTiming( time_index );
00198         hypre_PrintTiming( "Solve phase times", pcomm->comm() );
00199         hypre_FinalizeTiming( time_index );
00200         hypre_ClearTiming();
00201 #endif
00202         HYPRE_ParCSRPCGGetNumIterations( pcg_solver, &num_iterations );
00203         HYPRE_ParCSRPCGGetFinalRelativeResidualNorm( pcg_solver, &final_res_norm );
00204         MPI_Comm_rank( pcomm->comm(), &myid );
00205 
00206         if( myid == 0 )
00207         {
00208             cout << "PCG Iterations = " << num_iterations << endl
00209                  << "Final PCG Relative Residual Norm = " << final_res_norm << endl;
00210         }
00211     }
00212 
00213     return err;
00214 }
00215 
00216 HyprePCG::~HyprePCG()
00217 {
00218     HYPRE_ParCSRPCGDestroy( pcg_solver );
00219 }
00220 
00221 HypreGMRES::HypreGMRES( HypreParMatrix& _A ) : HypreSolver( &_A, true )
00222 {
00223     MPI_Comm comm;
00224     int k_dim      = 50;
00225     int max_iter   = 100;
00226     double tol     = 1e-6;
00227     iterative_mode = true;
00228     HYPRE_ParCSRMatrixGetComm( *A, &comm );
00229     HYPRE_ParCSRGMRESCreate( comm, &gmres_solver );
00230     HYPRE_ParCSRGMRESSetKDim( gmres_solver, k_dim );
00231     HYPRE_ParCSRGMRESSetMaxIter( gmres_solver, max_iter );
00232     HYPRE_ParCSRGMRESSetTol( gmres_solver, tol );
00233 }
00234 
00235 void HypreGMRES::SetTol( double atol, double rtol )
00236 {
00237     HYPRE_GMRESSetAbsoluteTol( gmres_solver, atol );
00238     HYPRE_GMRESSetTol( gmres_solver, rtol );
00239 }
00240 
00241 void HypreGMRES::SetMaxIter( int max_iter )
00242 {
00243     HYPRE_GMRESSetMaxIter( gmres_solver, max_iter );
00244 }
00245 
00246 void HypreGMRES::SetKDim( int k_dim )
00247 {
00248     HYPRE_GMRESSetKDim( gmres_solver, k_dim );
00249 }
00250 
00251 void HypreGMRES::SetLogging( int logging )
00252 {
00253     HYPRE_GMRESSetLogging( gmres_solver, logging );
00254 }
00255 
00256 void HypreGMRES::Verbosity( int print_lvl )
00257 {
00258     HYPRE_GMRESSetPrintLevel( gmres_solver, print_lvl );
00259 }
00260 
00261 void HypreGMRES::SetPreconditioner( HypreSolver& precond )
00262 {
00263     HYPRE_ParCSRGMRESSetPrecond( gmres_solver, precond.SolveFcn(), precond.SetupFcn(), precond );
00264 }
00265 
00266 HYPRE_Int HypreGMRES::Solve( const HypreParVector& b, HypreParVector& x ) const
00267 {
00268     HYPRE_Int err;
00269     int myid;
00270     HYPRE_Int time_index = 0;
00271     HYPRE_Int num_iterations;
00272     double final_res_norm;
00273     MPI_Comm comm;
00274     HYPRE_Int print_level;
00275     HYPRE_GMRESGetPrintLevel( gmres_solver, &print_level );
00276     HYPRE_ParCSRMatrixGetComm( *A, &comm );
00277 
00278     if( !setup_called )
00279     {
00280         if( print_level > 0 )
00281         {
00282 #ifdef HYPRE_TIMING
00283             time_index = hypre_InitializeTiming( "GMRES Setup" );
00284             hypre_BeginTiming( time_index );
00285 #endif
00286         }
00287 
00288         err          = HYPRE_ParCSRGMRESSetup( gmres_solver, *A, b, x );
00289         setup_called = 1;
00290 
00291         if( err != 0 )
00292         {
00293             cout << "PCG Setup failed with error = " << err << endl;
00294             return err;
00295         }
00296 
00297         if( print_level > 0 )
00298         {
00299 #ifdef HYPRE_TIMING
00300             hypre_EndTiming( time_index );
00301             hypre_PrintTiming( "Setup phase times", comm );
00302             hypre_FinalizeTiming( time_index );
00303             hypre_ClearTiming();
00304 #endif
00305         }
00306     }
00307 
00308     if( print_level > 0 )
00309     {
00310 #ifdef HYPRE_TIMING
00311         time_index = hypre_InitializeTiming( "GMRES Solve" );
00312         hypre_BeginTiming( time_index );
00313 #endif
00314     }
00315 
00316     if( !iterative_mode )
00317     {
00318         x = 0.0;
00319     }
00320 
00321     err = HYPRE_ParCSRGMRESSolve( gmres_solver, *A, b, x );
00322 
00323     if( print_level > 0 )
00324     {
00325 #ifdef HYPRE_TIMING
00326         hypre_EndTiming( time_index );
00327         hypre_PrintTiming( "Solve phase times", comm );
00328         hypre_FinalizeTiming( time_index );
00329         hypre_ClearTiming();
00330 #endif
00331         HYPRE_ParCSRGMRESGetNumIterations( gmres_solver, &num_iterations );
00332         HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm( gmres_solver, &final_res_norm );
00333         MPI_Comm_rank( comm, &myid );
00334 
00335         if( myid == 0 )
00336         {
00337             cout << "GMRES Iterations = " << num_iterations << endl
00338                  << "Final GMRES Relative Residual Norm = " << final_res_norm << endl;
00339         }
00340     }
00341 
00342     return err;
00343 }
00344 
00345 HypreGMRES::~HypreGMRES()
00346 {
00347     HYPRE_ParCSRGMRESDestroy( gmres_solver );
00348 }
00349 
00350 HypreParaSails::HypreParaSails( HypreParMatrix& A ) : HypreSolver( &A )
00351 {
00352     MPI_Comm comm;
00353     int sai_max_levels   = 1;
00354     double sai_threshold = 0.1;
00355     double sai_filter    = 0.1;
00356     int sai_sym          = 0;
00357     double sai_loadbal   = 0.0;
00358     int sai_reuse        = 0;
00359     int sai_logging      = 1;
00360     HYPRE_ParCSRMatrixGetComm( A, &comm );
00361     HYPRE_ParaSailsCreate( comm, &sai_precond );
00362     HYPRE_ParaSailsSetParams( sai_precond, sai_threshold, sai_max_levels );
00363     HYPRE_ParaSailsSetFilter( sai_precond, sai_filter );
00364     HYPRE_ParaSailsSetSym( sai_precond, sai_sym );
00365     HYPRE_ParaSailsSetLoadbal( sai_precond, sai_loadbal );
00366     HYPRE_ParaSailsSetReuse( sai_precond, sai_reuse );
00367     HYPRE_ParaSailsSetLogging( sai_precond, sai_logging );
00368 }
00369 
00370 void HypreParaSails::SetSymmetry( int sym )
00371 {
00372     HYPRE_ParaSailsSetSym( sai_precond, sym );
00373 }
00374 
00375 HypreParaSails::~HypreParaSails()
00376 {
00377     HYPRE_ParaSailsDestroy( sai_precond );
00378 }
00379 
00380 HypreBoomerAMG::HypreBoomerAMG()
00381 {
00382     HYPRE_BoomerAMGCreate( &amg_precond );
00383     SetDefaultOptions();
00384 }
00385 
00386 HypreBoomerAMG::HypreBoomerAMG( HypreParMatrix& A ) : HypreSolver( &A )
00387 {
00388     HYPRE_BoomerAMGCreate( &amg_precond );
00389     SetDefaultOptions();
00390 }
00391 
00392 void HypreBoomerAMG::SetDefaultOptions()
00393 {
00394     // AMG coarsening options:
00395     int coarsen_type = 10;    // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
00396     int agg_levels   = 2;     // number of aggressive coarsening levels
00397     double theta     = 0.25;  // strength threshold: 0.25, 0.5, 0.8
00398     // AMG interpolation options:
00399     int interp_type = 6;  // 6 = extended+i, 0 = classical
00400     int Pmax        = 4;  // max number of elements per row in P
00401     // AMG relaxation options:
00402     int relax_type   = 10;  // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
00403     int relax_sweeps = 1;   // relaxation sweeps on each level
00404     // Additional options:
00405     int print_level = 1;   // print AMG iterations? 1 = no, 2 = yes
00406     int max_levels  = 25;  // max number of levels in AMG hierarchy
00407     HYPRE_BoomerAMGSetCoarsenType( amg_precond, coarsen_type );
00408     HYPRE_BoomerAMGSetAggNumLevels( amg_precond, agg_levels );
00409     HYPRE_BoomerAMGSetRelaxType( amg_precond, relax_type );
00410     HYPRE_BoomerAMGSetNumSweeps( amg_precond, relax_sweeps );
00411     HYPRE_BoomerAMGSetStrongThreshold( amg_precond, theta );
00412     HYPRE_BoomerAMGSetInterpType( amg_precond, interp_type );
00413     HYPRE_BoomerAMGSetPMaxElmts( amg_precond, Pmax );
00414     HYPRE_BoomerAMGSetPrintLevel( amg_precond, print_level );
00415     HYPRE_BoomerAMGSetMaxLevels( amg_precond, max_levels );
00416     // Use as a preconditioner (one V-cycle, zero tolerance)
00417     HYPRE_BoomerAMGSetMaxIter( amg_precond, 1 );
00418     HYPRE_BoomerAMGSetTol( amg_precond, 0.0 );
00419 }
00420 
00421 void HypreBoomerAMG::ResetAMGPrecond()
00422 {
00423     HYPRE_Int coarsen_type;
00424     HYPRE_Int agg_levels;
00425     HYPRE_Int relax_type;
00426     HYPRE_Int relax_sweeps;
00427     HYPRE_Real theta;
00428     HYPRE_Int interp_type;
00429     HYPRE_Int Pmax;
00430     HYPRE_Int print_level;
00431     HYPRE_Int dim;
00432     hypre_ParAMGData* amg_data = (hypre_ParAMGData*)amg_precond;
00433     // read options from amg_precond
00434     HYPRE_BoomerAMGGetCoarsenType( amg_precond, &coarsen_type );
00435     agg_levels   = hypre_ParAMGDataAggNumLevels( amg_data );
00436     relax_type   = hypre_ParAMGDataUserRelaxType( amg_data );
00437     relax_sweeps = hypre_ParAMGDataUserNumSweeps( amg_data );
00438     HYPRE_BoomerAMGGetStrongThreshold( amg_precond, &theta );
00439     hypre_BoomerAMGGetInterpType( amg_precond, &interp_type );
00440     HYPRE_BoomerAMGGetPMaxElmts( amg_precond, &Pmax );
00441     HYPRE_BoomerAMGGetPrintLevel( amg_precond, &print_level );
00442     HYPRE_BoomerAMGGetNumFunctions( amg_precond, &dim );
00443     HYPRE_BoomerAMGDestroy( amg_precond );
00444     HYPRE_BoomerAMGCreate( &amg_precond );
00445     HYPRE_BoomerAMGSetCoarsenType( amg_precond, coarsen_type );
00446     HYPRE_BoomerAMGSetAggNumLevels( amg_precond, agg_levels );
00447     HYPRE_BoomerAMGSetRelaxType( amg_precond, relax_type );
00448     HYPRE_BoomerAMGSetNumSweeps( amg_precond, relax_sweeps );
00449     HYPRE_BoomerAMGSetMaxLevels( amg_precond, 25 );
00450     HYPRE_BoomerAMGSetTol( amg_precond, 0.0 );
00451     HYPRE_BoomerAMGSetMaxIter( amg_precond, 1 );  // one V-cycle
00452     HYPRE_BoomerAMGSetStrongThreshold( amg_precond, theta );
00453     HYPRE_BoomerAMGSetInterpType( amg_precond, interp_type );
00454     HYPRE_BoomerAMGSetPMaxElmts( amg_precond, Pmax );
00455     HYPRE_BoomerAMGSetPrintLevel( amg_precond, print_level );
00456     HYPRE_BoomerAMGSetNumFunctions( amg_precond, dim );
00457 }
00458 
00459 void HypreBoomerAMG::SetOperator( const HypreParMatrix& op )
00460 {
00461     if( A )
00462     {
00463         ResetAMGPrecond();
00464     }
00465 
00466     // update base classes: Operator, Solver, HypreSolver
00467     A            = const_cast< HypreParMatrix* >( &op );
00468     setup_called = 0;
00469     delete X;
00470     delete B;
00471     B = X = NULL;
00472 }
00473 
00474 void HypreBoomerAMG::SetSystemsOptions( int dim )
00475 {
00476     HYPRE_BoomerAMGSetNumFunctions( amg_precond, dim );
00477     // More robust options with respect to convergence
00478     HYPRE_BoomerAMGSetAggNumLevels( amg_precond, 0 );
00479     HYPRE_BoomerAMGSetStrongThreshold( amg_precond, dim * 0.25 );
00480 }
00481 
00482 HypreBoomerAMG::~HypreBoomerAMG()
00483 {
00484     HYPRE_BoomerAMGDestroy( amg_precond );
00485 }
00486 
00487 }  // namespace moab
00488 
00489 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines