Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
HypreSolver.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "HypreSolver.hpp"
13 #include "HypreParMatrix.hpp"
14 
15 #ifdef MOAB_HAVE_MPI
16 
17 #include "hypre_parcsr.hpp"
18 
19 #include <fstream>
20 #include <iostream>
21 #include <cmath>
22 #include <cstdlib>
23 #include <cassert>
24 
25 using namespace std;
26 
27 namespace moab
28 {
29 #define moab_hypre_assert( a, b ) \
30  { \
31  }
32 #define moab_hypre_assert_t( a, b ) \
33  { \
34  if( !a ) \
35  { \
36  std::cout << "HYPRE Error: " << b << std::endl; \
37  exit( -1 ); \
38  }
39 
40 HypreSolver::HypreSolver( bool iterative ) : AbstractSolver( iterative )
41 {
42  A = NULL;
43  setup_called = 0;
44  B = X = NULL;
45 }
46 
47 HypreSolver::HypreSolver( HypreParMatrix* _A, bool iterative ) : AbstractSolver( iterative )
48 {
49  A = _A;
50  setup_called = 0;
51  B = X = NULL;
52 }
53 
54 HYPRE_Int HypreSolver::Solve( const HypreParVector& b, HypreParVector& x ) const
55 {
56  HYPRE_Int err;
57 
58  if( A == NULL )
59  {
60  MB_SET_ERR_RET_VAL( "HypreSolver::Solve (...) : HypreParMatrix A is missing", 1 );
61  }
62 
63  if( !setup_called )
64  {
65  err = SetupFcn()( *this, *A, b.x_par, x.x_par );
66  setup_called = 1;
67  }
68 
69  if( !iterative_mode )
70  {
71  x = 0.0;
72  }
73 
74  err = SolveFcn()( *this, *A, b.x_par, x.x_par );
75  return err;
76 }
77 
78 HypreSolver::~HypreSolver()
79 {
80  if( B )
81  {
82  delete B;
83  }
84 
85  if( X )
86  {
87  delete X;
88  }
89 }
90 
91 HyprePCG::HyprePCG( HypreParMatrix& _A ) : HypreSolver( &_A, true )
92 {
93  iterative_mode = true;
94  HYPRE_ParCSRPCGCreate( A->GetParallelCommunicator()->comm(), &pcg_solver );
95 }
96 
97 void HyprePCG::SetTol( double tol )
98 {
99  HYPRE_PCGSetTol( pcg_solver, tol );
100 }
101 
102 void HyprePCG::SetMaxIter( int max_iter )
103 {
104  HYPRE_PCGSetMaxIter( pcg_solver, max_iter );
105 }
106 
107 void HyprePCG::SetLogging( int logging )
108 {
109  HYPRE_PCGSetLogging( pcg_solver, logging );
110 }
111 
112 void HyprePCG::Verbosity( int print_lvl )
113 {
114  HYPRE_ParCSRPCGSetPrintLevel( pcg_solver, print_lvl );
115 }
116 
117 void HyprePCG::SetPreconditioner( HypreSolver& precond )
118 {
119  HYPRE_ParCSRPCGSetPrecond( pcg_solver, precond.SolveFcn(), precond.SetupFcn(), precond );
120 }
121 
122 void HyprePCG::SetResidualConvergenceOptions( int res_frequency, double rtol )
123 {
124  HYPRE_PCGSetTwoNorm( pcg_solver, 1 );
125 
126  if( res_frequency > 0 )
127  {
128  HYPRE_PCGSetRecomputeResidualP( pcg_solver, res_frequency );
129  }
130 
131  if( rtol > 0.0 )
132  {
133  HYPRE_PCGSetResidualTol( pcg_solver, rtol );
134  }
135 }
136 
137 HYPRE_Int HyprePCG::Solve( const HypreParVector& b, HypreParVector& x ) const
138 {
139  HYPRE_Int err;
140  int myid;
141  HYPRE_Int time_index = 0;
142  HYPRE_Int num_iterations;
143  double final_res_norm;
144  moab::ParallelComm* pcomm;
145  HYPRE_Int print_level;
146  HYPRE_PCGGetPrintLevel( pcg_solver, &print_level );
147  pcomm = A->GetParallelCommunicator();
148 
149  if( !setup_called )
150  {
151  if( print_level > 0 )
152  {
153 #ifdef HYPRE_TIMING
154  time_index = hypre_InitializeTiming( "PCG Setup" );
155  hypre_BeginTiming( time_index );
156 #endif
157  }
158 
159  err = HYPRE_ParCSRPCGSetup( pcg_solver, *A, (HYPRE_ParVector)b, (HYPRE_ParVector)x );
160  setup_called = 1;
161 
162  if( err != 0 )
163  {
164  cout << "PCG Setup failed with error = " << err << endl;
165  return err;
166  }
167 
168  if( print_level > 0 )
169  {
170 #ifdef HYPRE_TIMING
171  hypre_EndTiming( time_index );
172  hypre_PrintTiming( "Setup phase times", pcomm->comm() );
173  hypre_FinalizeTiming( time_index );
174  hypre_ClearTiming();
175 #endif
176  }
177  }
178 
179  if( print_level > 0 )
180  {
181 #ifdef HYPRE_TIMING
182  time_index = hypre_InitializeTiming( "PCG Solve" );
183  hypre_BeginTiming( time_index );
184 #endif
185  }
186 
187  if( !iterative_mode )
188  {
189  x = 0.0;
190  }
191 
192  err = HYPRE_ParCSRPCGSolve( pcg_solver, *A, (HYPRE_ParVector)b, (HYPRE_ParVector)x );
193 
194  if( print_level > 0 )
195  {
196 #ifdef HYPRE_TIMING
197  hypre_EndTiming( time_index );
198  hypre_PrintTiming( "Solve phase times", pcomm->comm() );
199  hypre_FinalizeTiming( time_index );
200  hypre_ClearTiming();
201 #endif
202  HYPRE_ParCSRPCGGetNumIterations( pcg_solver, &num_iterations );
203  HYPRE_ParCSRPCGGetFinalRelativeResidualNorm( pcg_solver, &final_res_norm );
204  MPI_Comm_rank( pcomm->comm(), &myid );
205 
206  if( myid == 0 )
207  {
208  cout << "PCG Iterations = " << num_iterations << endl
209  << "Final PCG Relative Residual Norm = " << final_res_norm << endl;
210  }
211  }
212 
213  return err;
214 }
215 
216 HyprePCG::~HyprePCG()
217 {
218  HYPRE_ParCSRPCGDestroy( pcg_solver );
219 }
220 
221 HypreGMRES::HypreGMRES( HypreParMatrix& _A ) : HypreSolver( &_A, true )
222 {
223  MPI_Comm comm;
224  int k_dim = 50;
225  int max_iter = 100;
226  double tol = 1e-6;
227  iterative_mode = true;
228  HYPRE_ParCSRMatrixGetComm( *A, &comm );
229  HYPRE_ParCSRGMRESCreate( comm, &gmres_solver );
230  HYPRE_ParCSRGMRESSetKDim( gmres_solver, k_dim );
231  HYPRE_ParCSRGMRESSetMaxIter( gmres_solver, max_iter );
232  HYPRE_ParCSRGMRESSetTol( gmres_solver, tol );
233 }
234 
235 void HypreGMRES::SetTol( double atol, double rtol )
236 {
237  HYPRE_GMRESSetAbsoluteTol( gmres_solver, atol );
238  HYPRE_GMRESSetTol( gmres_solver, rtol );
239 }
240 
241 void HypreGMRES::SetMaxIter( int max_iter )
242 {
243  HYPRE_GMRESSetMaxIter( gmres_solver, max_iter );
244 }
245 
246 void HypreGMRES::SetKDim( int k_dim )
247 {
248  HYPRE_GMRESSetKDim( gmres_solver, k_dim );
249 }
250 
251 void HypreGMRES::SetLogging( int logging )
252 {
253  HYPRE_GMRESSetLogging( gmres_solver, logging );
254 }
255 
256 void HypreGMRES::Verbosity( int print_lvl )
257 {
258  HYPRE_GMRESSetPrintLevel( gmres_solver, print_lvl );
259 }
260 
261 void HypreGMRES::SetPreconditioner( HypreSolver& precond )
262 {
263  HYPRE_ParCSRGMRESSetPrecond( gmres_solver, precond.SolveFcn(), precond.SetupFcn(), precond );
264 }
265 
266 HYPRE_Int HypreGMRES::Solve( const HypreParVector& b, HypreParVector& x ) const
267 {
268  HYPRE_Int err;
269  int myid;
270  HYPRE_Int time_index = 0;
271  HYPRE_Int num_iterations;
272  double final_res_norm;
273  MPI_Comm comm;
274  HYPRE_Int print_level;
275  HYPRE_GMRESGetPrintLevel( gmres_solver, &print_level );
276  HYPRE_ParCSRMatrixGetComm( *A, &comm );
277 
278  if( !setup_called )
279  {
280  if( print_level > 0 )
281  {
282 #ifdef HYPRE_TIMING
283  time_index = hypre_InitializeTiming( "GMRES Setup" );
284  hypre_BeginTiming( time_index );
285 #endif
286  }
287 
288  err = HYPRE_ParCSRGMRESSetup( gmres_solver, *A, b, x );
289  setup_called = 1;
290 
291  if( err != 0 )
292  {
293  cout << "PCG Setup failed with error = " << err << endl;
294  return err;
295  }
296 
297  if( print_level > 0 )
298  {
299 #ifdef HYPRE_TIMING
300  hypre_EndTiming( time_index );
301  hypre_PrintTiming( "Setup phase times", comm );
302  hypre_FinalizeTiming( time_index );
303  hypre_ClearTiming();
304 #endif
305  }
306  }
307 
308  if( print_level > 0 )
309  {
310 #ifdef HYPRE_TIMING
311  time_index = hypre_InitializeTiming( "GMRES Solve" );
312  hypre_BeginTiming( time_index );
313 #endif
314  }
315 
316  if( !iterative_mode )
317  {
318  x = 0.0;
319  }
320 
321  err = HYPRE_ParCSRGMRESSolve( gmres_solver, *A, b, x );
322 
323  if( print_level > 0 )
324  {
325 #ifdef HYPRE_TIMING
326  hypre_EndTiming( time_index );
327  hypre_PrintTiming( "Solve phase times", comm );
328  hypre_FinalizeTiming( time_index );
329  hypre_ClearTiming();
330 #endif
331  HYPRE_ParCSRGMRESGetNumIterations( gmres_solver, &num_iterations );
332  HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm( gmres_solver, &final_res_norm );
333  MPI_Comm_rank( comm, &myid );
334 
335  if( myid == 0 )
336  {
337  cout << "GMRES Iterations = " << num_iterations << endl
338  << "Final GMRES Relative Residual Norm = " << final_res_norm << endl;
339  }
340  }
341 
342  return err;
343 }
344 
345 HypreGMRES::~HypreGMRES()
346 {
347  HYPRE_ParCSRGMRESDestroy( gmres_solver );
348 }
349 
350 HypreParaSails::HypreParaSails( HypreParMatrix& A ) : HypreSolver( &A )
351 {
352  MPI_Comm comm;
353  int sai_max_levels = 1;
354  double sai_threshold = 0.1;
355  double sai_filter = 0.1;
356  int sai_sym = 0;
357  double sai_loadbal = 0.0;
358  int sai_reuse = 0;
359  int sai_logging = 1;
360  HYPRE_ParCSRMatrixGetComm( A, &comm );
361  HYPRE_ParaSailsCreate( comm, &sai_precond );
362  HYPRE_ParaSailsSetParams( sai_precond, sai_threshold, sai_max_levels );
363  HYPRE_ParaSailsSetFilter( sai_precond, sai_filter );
364  HYPRE_ParaSailsSetSym( sai_precond, sai_sym );
365  HYPRE_ParaSailsSetLoadbal( sai_precond, sai_loadbal );
366  HYPRE_ParaSailsSetReuse( sai_precond, sai_reuse );
367  HYPRE_ParaSailsSetLogging( sai_precond, sai_logging );
368 }
369 
370 void HypreParaSails::SetSymmetry( int sym )
371 {
372  HYPRE_ParaSailsSetSym( sai_precond, sym );
373 }
374 
375 HypreParaSails::~HypreParaSails()
376 {
377  HYPRE_ParaSailsDestroy( sai_precond );
378 }
379 
380 HypreBoomerAMG::HypreBoomerAMG()
381 {
382  HYPRE_BoomerAMGCreate( &amg_precond );
383  SetDefaultOptions();
384 }
385 
386 HypreBoomerAMG::HypreBoomerAMG( HypreParMatrix& A ) : HypreSolver( &A )
387 {
388  HYPRE_BoomerAMGCreate( &amg_precond );
389  SetDefaultOptions();
390 }
391 
392 void HypreBoomerAMG::SetDefaultOptions()
393 {
394  // AMG coarsening options:
395  int coarsen_type = 10; // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
396  int agg_levels = 2; // number of aggressive coarsening levels
397  double theta = 0.25; // strength threshold: 0.25, 0.5, 0.8
398  // AMG interpolation options:
399  int interp_type = 6; // 6 = extended+i, 0 = classical
400  int Pmax = 4; // max number of elements per row in P
401  // AMG relaxation options:
402  int relax_type = 10; // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
403  int relax_sweeps = 1; // relaxation sweeps on each level
404  // Additional options:
405  int print_level = 1; // print AMG iterations? 1 = no, 2 = yes
406  int max_levels = 25; // max number of levels in AMG hierarchy
407  HYPRE_BoomerAMGSetCoarsenType( amg_precond, coarsen_type );
408  HYPRE_BoomerAMGSetAggNumLevels( amg_precond, agg_levels );
409  HYPRE_BoomerAMGSetRelaxType( amg_precond, relax_type );
410  HYPRE_BoomerAMGSetNumSweeps( amg_precond, relax_sweeps );
411  HYPRE_BoomerAMGSetStrongThreshold( amg_precond, theta );
412  HYPRE_BoomerAMGSetInterpType( amg_precond, interp_type );
413  HYPRE_BoomerAMGSetPMaxElmts( amg_precond, Pmax );
414  HYPRE_BoomerAMGSetPrintLevel( amg_precond, print_level );
415  HYPRE_BoomerAMGSetMaxLevels( amg_precond, max_levels );
416  // Use as a preconditioner (one V-cycle, zero tolerance)
417  HYPRE_BoomerAMGSetMaxIter( amg_precond, 1 );
418  HYPRE_BoomerAMGSetTol( amg_precond, 0.0 );
419 }
420 
421 void HypreBoomerAMG::ResetAMGPrecond()
422 {
423  HYPRE_Int coarsen_type;
424  HYPRE_Int agg_levels;
425  HYPRE_Int relax_type;
426  HYPRE_Int relax_sweeps;
427  HYPRE_Real theta;
428  HYPRE_Int interp_type;
429  HYPRE_Int Pmax;
430  HYPRE_Int print_level;
431  HYPRE_Int dim;
432  hypre_ParAMGData* amg_data = (hypre_ParAMGData*)amg_precond;
433  // read options from amg_precond
434  HYPRE_BoomerAMGGetCoarsenType( amg_precond, &coarsen_type );
435  agg_levels = hypre_ParAMGDataAggNumLevels( amg_data );
436  relax_type = hypre_ParAMGDataUserRelaxType( amg_data );
437  relax_sweeps = hypre_ParAMGDataUserNumSweeps( amg_data );
438  HYPRE_BoomerAMGGetStrongThreshold( amg_precond, &theta );
439  hypre_BoomerAMGGetInterpType( amg_precond, &interp_type );
440  HYPRE_BoomerAMGGetPMaxElmts( amg_precond, &Pmax );
441  HYPRE_BoomerAMGGetPrintLevel( amg_precond, &print_level );
442  HYPRE_BoomerAMGGetNumFunctions( amg_precond, &dim );
443  HYPRE_BoomerAMGDestroy( amg_precond );
444  HYPRE_BoomerAMGCreate( &amg_precond );
445  HYPRE_BoomerAMGSetCoarsenType( amg_precond, coarsen_type );
446  HYPRE_BoomerAMGSetAggNumLevels( amg_precond, agg_levels );
447  HYPRE_BoomerAMGSetRelaxType( amg_precond, relax_type );
448  HYPRE_BoomerAMGSetNumSweeps( amg_precond, relax_sweeps );
449  HYPRE_BoomerAMGSetMaxLevels( amg_precond, 25 );
450  HYPRE_BoomerAMGSetTol( amg_precond, 0.0 );
451  HYPRE_BoomerAMGSetMaxIter( amg_precond, 1 ); // one V-cycle
452  HYPRE_BoomerAMGSetStrongThreshold( amg_precond, theta );
453  HYPRE_BoomerAMGSetInterpType( amg_precond, interp_type );
454  HYPRE_BoomerAMGSetPMaxElmts( amg_precond, Pmax );
455  HYPRE_BoomerAMGSetPrintLevel( amg_precond, print_level );
456  HYPRE_BoomerAMGSetNumFunctions( amg_precond, dim );
457 }
458 
459 void HypreBoomerAMG::SetOperator( const HypreParMatrix& op )
460 {
461  if( A )
462  {
463  ResetAMGPrecond();
464  }
465 
466  // update base classes: Operator, Solver, HypreSolver
467  A = const_cast< HypreParMatrix* >( &op );
468  setup_called = 0;
469  delete X;
470  delete B;
471  B = X = NULL;
472 }
473 
474 void HypreBoomerAMG::SetSystemsOptions( int dim )
475 {
476  HYPRE_BoomerAMGSetNumFunctions( amg_precond, dim );
477  // More robust options with respect to convergence
478  HYPRE_BoomerAMGSetAggNumLevels( amg_precond, 0 );
479  HYPRE_BoomerAMGSetStrongThreshold( amg_precond, dim * 0.25 );
480 }
481 
482 HypreBoomerAMG::~HypreBoomerAMG()
483 {
484  HYPRE_BoomerAMGDestroy( amg_precond );
485 }
486 
487 } // namespace moab
488 
489 #endif