 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
1.9.1.