hypre_solver.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1132 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-02-23 05:43:26 +0000 (Tue, 23 Feb 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #include "hypre_solver.h"
31 
32 // For problem->get_jacobian(...)
33 #include "problem.h"
34 
35 
36 namespace oomph
37 {
38 
39 ////////////////////////////////////////////////////////////////////////
40 ////////////////////////////////////////////////////////////////////////
41 ////////////////////////////////////////////////////////////////////////
42 
43 //==================================================================
44 /// Default settings for various uses of the HYPRE solver
45 //==================================================================
46  namespace Hypre_default_settings
47  {
48  /// \short Set default parameters for use as preconditioner in
49  /// 2D Poisson-type problem.
51  HyprePreconditioner* hypre_preconditioner_pt)
52  {
53  // Set iterations to 1
54  hypre_preconditioner_pt->set_amg_iterations(1);
55 
56  // Use simple smoother
57  hypre_preconditioner_pt->amg_using_simple_smoothing();
58 
59  // Smoother types:
60  // 0=Jacobi
61  // 1=Gauss-Seidel
62  hypre_preconditioner_pt->amg_simple_smoother() = 1;
63 
64  // AMG preconditioner
65  hypre_preconditioner_pt->hypre_method() = HyprePreconditioner::BoomerAMG;
66 
67  // Choose strength parameter for amg
68  hypre_preconditioner_pt->amg_strength() = 0.25;
69 
70  // Coarsening type
71  hypre_preconditioner_pt->amg_coarsening() = 0;
72  }
73 
74 
75  /// \short Set default parameters for use as preconditioner in
76  /// 3D Poisson-type problem.
78  HyprePreconditioner* hypre_preconditioner_pt)
79  {
80  // Set default settings as for 2D Poisson
81  set_defaults_for_2D_poisson_problem(hypre_preconditioner_pt);
82 
83  // Change strength parameter for amg
84  hypre_preconditioner_pt->amg_strength() = 0.7;
85  }
86 
87 
88  /// \short Set default parameters for use as preconditioner in
89  /// for momentum block in Navier-Stokes problem
91  HyprePreconditioner* hypre_preconditioner_pt)
92  {
93  // Set default settings as for 2D Poisson
94  set_defaults_for_2D_poisson_problem(hypre_preconditioner_pt);
95 
96  // Change smoother type:
97  // 0=Jacobi
98  // 1=Gauss-Seidel
99  hypre_preconditioner_pt->amg_simple_smoother() = 0;
100 
101  // Set smoother damping
102  hypre_preconditioner_pt->amg_damping() = 0.5;
103 
104  // Change strength parameter for amg
105  hypre_preconditioner_pt->amg_strength() = 0.75;
106  }
107 
108  }
109 
110 
111 ////////////////////////////////////////////////////////////////////////
112 ////////////////////////////////////////////////////////////////////////
113 // functions for HypreHelpers namespace
114 ////////////////////////////////////////////////////////////////////////
115 ////////////////////////////////////////////////////////////////////////
116 
117  namespace HypreHelpers
118  {
119 
120 
121 //========================================================================
122  /// \short Default for AMG strength (0.25 recommended for 2D problems;
123  /// larger (0.5-0.75, say) for 3D
124 //========================================================================
125  double AMG_strength=0.25;
126 
127 //========================================================================
128  /// \short Default AMG coarsening strategy. Coarsening types include:
129  /// 0 = CLJP (parallel coarsening using independent sets)
130  /// 1 = classical RS with no boundary treatment (not recommended
131  /// in parallel)
132  /// 3 = modified RS with 3rd pass to add C points on the boundaries
133  /// 6 = Falgout (uses 1 then CLJP using interior coarse points as
134  /// first independent set)
135  /// 8 = PMIS (parallel coarsening using independent sets - lower
136  /// complexities than 0, maybe also slower convergence)
137  /// 10= HMIS (one pass RS on each processor then PMIS on interior
138  /// coarse points as first independent set)
139  /// 11= One pass RS on each processor (not recommended)
140 //========================================================================
141  unsigned AMG_coarsening=6;
142 
143 
144 //========================================================================
145  /// AMG interpolation truncation factor
146 //========================================================================
147  double AMG_truncation=0.0;
148 
149 //========================================================================
150 /// Helper function to check the Hypre error flag, return the message
151 /// associated with any error, and reset the global error flag to zero.
152 /// This function also returns the error value.
153 //========================================================================
154  int check_HYPRE_error_flag(std::ostringstream& message)
155  {
156  // get the Hypre error flag
157  int err = HYPRE_GetError();
158 
159  // Tell us all about it...
160  if (err)
161  {
162  oomph_info << "Hypre error flag=" << err << std::endl;
163  char* error_message = new char[128];
164  HYPRE_DescribeError(err, error_message);
165  message << "WARNING: " << std::endl
166  << "HYPRE error message: " << error_message
167  << std::endl;
168  delete[] error_message;
169  }
170 
171  // reset Hypre's global error flag
173 
174  return err;
175  }
176 
177 
178 //========================================================================
179 /// Helper function to create a HYPRE_IJVector and HYPRE_ParVector.
180 /// + If no MPI then serial vectors are created
181 /// + If MPI and serial input vector then distributed hypre vectors are created
182 /// + If MPI and distributed input vector the distributed output vectors
183 /// are created.
184 //========================================================================
185  void create_HYPRE_Vector(const DoubleVector &oomph_vec, const
186  LinearAlgebraDistribution* dist_pt,
187  HYPRE_IJVector& hypre_ij_vector,
188  HYPRE_ParVector& hypre_par_vector)
189  {
190  // the lower and upper row of the vector on this processor
191  unsigned lower = dist_pt->first_row();
192  unsigned upper = dist_pt->first_row() +
193  dist_pt->nrow_local() - 1;
194 
195  // number of local rows
196  unsigned nrow_local = dist_pt->nrow_local();
197 
198  // initialize Hypre vector
199 #ifdef OOMPH_HAS_MPI
200  HYPRE_IJVectorCreate(oomph_vec.distribution_pt()->
201  communicator_pt()->mpi_comm(),
202  lower,upper, &hypre_ij_vector);
203 #else
204  HYPRE_IJVectorCreate(MPI_COMM_WORLD,
205  lower,upper, &hypre_ij_vector);
206 #endif
207  HYPRE_IJVectorSetObjectType(hypre_ij_vector, HYPRE_PARCSR);
208  HYPRE_IJVectorInitialize(hypre_ij_vector);
209 
210  // set up array containing indices
211  int* indices = new int[nrow_local];
212  double* values = new double[nrow_local];
213  const unsigned hypre_first_row = dist_pt->first_row();
214  unsigned j = 0;
215  if (!oomph_vec.distributed() && dist_pt->distributed())
216  {
217  j = hypre_first_row;
218  }
219  const double* o_pt = oomph_vec.values_pt();
220  for (unsigned i=0; i<nrow_local; i++)
221  {
222  indices[i] = hypre_first_row + i;
223  values[i] = o_pt[j+i];
224  }
225 
226  // insert values
227  HYPRE_IJVectorSetValues(hypre_ij_vector,nrow_local,indices,values);
228 
229  // assemble vectors
230  HYPRE_IJVectorAssemble(hypre_ij_vector);
231  HYPRE_IJVectorGetObject(hypre_ij_vector, (void **) &hypre_par_vector);
232 
233  // clean up
234  delete[] indices;
235  delete[] values;
236  }
237 
238 
239 //========================================================================
240 /// Helper function to create a HYPRE_IJVector and HYPRE_ParVector.
241 /// + If no MPI then serial vectors are created
242 /// + If MPI and serial input vector then distributed hypre vectors are created
243 /// + If MPI and distributed input vector the distributed output vectors
244 /// are created.
245 //========================================================================
247  HYPRE_IJVector& hypre_ij_vector,
248  HYPRE_ParVector& hypre_par_vector)
249  {
250  // the lower and upper row of the vector on this processor
251  unsigned lower = dist_pt->first_row();
252  unsigned upper = dist_pt->first_row() +
253  dist_pt->nrow_local() - 1;
254 
255  // initialize Hypre vector
256 #ifdef OOMPH_HAS_MPI
257  HYPRE_IJVectorCreate(dist_pt->communicator_pt()->mpi_comm(),
258  lower,upper, &hypre_ij_vector);
259 #else
260  HYPRE_IJVectorCreate(MPI_COMM_WORLD,
261  lower,upper, &hypre_ij_vector);
262 #endif
263  HYPRE_IJVectorSetObjectType(hypre_ij_vector, HYPRE_PARCSR);
264  HYPRE_IJVectorInitialize(hypre_ij_vector);
265 
266  // assemble vectors
267  HYPRE_IJVectorAssemble(hypre_ij_vector);
268  HYPRE_IJVectorGetObject(hypre_ij_vector, (void **) &hypre_par_vector);
269  }
270 
271 
272 //========================================================================
273 /// Helper function to create a serial HYPRE_IJMatrix and
274 /// HYPRE_ParCSRMatrix from a CRDoubleMatrix
275 /// NOTE: dist_pt is rebuilt to match the distribution of the hypre solver
276 /// which is not necassarily the same as the oomph lib matrix
277 //========================================================================
279  HYPRE_IJMatrix& hypre_ij_matrix,
280  HYPRE_ParCSRMatrix& hypre_par_matrix,
281  LinearAlgebraDistribution* dist_pt)
282  {
283 #ifdef PARANOID
284  // check that the matrix is built
285  if ( !oomph_matrix->built() )
286  {
287  std::ostringstream error_message;
288  error_message << "The matrix has not been built";
289  throw OomphLibError(error_message.str(),
290  OOMPH_CURRENT_FUNCTION,
291  OOMPH_EXCEPTION_LOCATION);
292  }
293  // check the matrix is square
294  if ( oomph_matrix->nrow() != oomph_matrix->ncol() )
295  {
296  std::ostringstream error_message;
297  error_message << "create_HYPRE_Matrix require a square matrix. "
298  << "Matrix is " << oomph_matrix->nrow()
299  << " by " << oomph_matrix->ncol() << std::endl;
300  throw OomphLibError(error_message.str(),
301  OOMPH_CURRENT_FUNCTION,
302  OOMPH_EXCEPTION_LOCATION);
303  }
304 #endif
305 
306 #ifdef OOMPH_HAS_MPI
307  //Trap the case when we have compiled with MPI,
308  //but are running in serial
310  {
311  std::ostringstream error_stream;
312  error_stream << "Oomph-lib has been compiled with MPI support and "
313  << "you are using HYPRE.\n"
314  <<
315  "For this combination of flags, MPI must be initialised.\n"
316  << "Call MPI_Helpers::init() in the "
317  << "main() function of your driver code\n";
318  throw OomphLibError(error_stream.str(),
319  OOMPH_CURRENT_FUNCTION,
320  OOMPH_EXCEPTION_LOCATION);
321  }
322 #endif
323 
324  // find number of rows/columns
325  const unsigned nrow = int(oomph_matrix->nrow());
326 
327  // get pointers to the matrix
328 
329  // column indices of matrix
330  const int* matrix_cols = oomph_matrix->column_index();
331 
332  // entries of matrix
333  const double* matrix_vals = oomph_matrix->value();
334 
335  // row starts
336  const int* matrix_row_start = oomph_matrix->row_start();
337 
338  // build the distribution
339  if (oomph_matrix->distribution_pt()->distributed())
340  {
341  dist_pt->build(oomph_matrix->distribution_pt());
342  }
343  else
344  {
345  bool distributed = true;
346  if (oomph_matrix->distribution_pt()->communicator_pt()->nproc() == 1)
347  {
348  distributed = false;
349  }
350  dist_pt->build
351  (oomph_matrix->distribution_pt()->communicator_pt(),nrow,distributed);
352  }
353 
354  // initialize hypre matrix
355  unsigned lower = dist_pt->first_row();
356  unsigned upper = lower + dist_pt->nrow_local() - 1;
357 
358 #ifdef OOMPH_HAS_MPI
359  HYPRE_IJMatrixCreate(dist_pt->communicator_pt()->mpi_comm(),
360  lower,
361  upper,
362  lower,
363  upper,
364  &hypre_ij_matrix);
365 #else
366  HYPRE_IJMatrixCreate(MPI_COMM_WORLD,
367  lower,
368  upper,
369  lower,
370  upper,
371  &hypre_ij_matrix);
372 #endif
373  HYPRE_IJMatrixSetObjectType(hypre_ij_matrix, HYPRE_PARCSR);
374  HYPRE_IJMatrixInitialize(hypre_ij_matrix);
375 
376  // set up a row map
377  // and first row / nrow_local
378  const unsigned hypre_nrow_local = dist_pt->nrow_local();
379  const unsigned hypre_first_row = dist_pt->first_row();
380  int* ncols_per_row = new int[hypre_nrow_local];
381  int* row_map = new int[hypre_nrow_local];
382  for (unsigned i = 0; i < hypre_nrow_local; i++)
383  {
384  unsigned j = i;
385  if (!oomph_matrix->distributed() && dist_pt->distributed())
386  {
387  j += hypre_first_row;
388  }
389  ncols_per_row[i] = matrix_row_start[j+1] - matrix_row_start[j];
390  row_map[i] = hypre_first_row + i;
391  }
392 
393  // put values in HYPRE matrix
394  int local_start = 0;
395  if (!oomph_matrix->distributed() && dist_pt->distributed())
396  {
397  local_start += matrix_row_start[hypre_first_row];
398  }
399 
400 
401  HYPRE_IJMatrixSetValues(hypre_ij_matrix,
402  hypre_nrow_local,
403  ncols_per_row,
404  row_map,
405  matrix_cols+local_start,
406  matrix_vals+local_start);
407 
408  // assemble matrix
409  HYPRE_IJMatrixAssemble(hypre_ij_matrix); // hierher leak?
410  HYPRE_IJMatrixGetObject(hypre_ij_matrix, (void **) &hypre_par_matrix);
411 
412  // tidy up memory
413  delete[] ncols_per_row;
414  delete[] row_map;
415  }
416 
417  //====================================================================
418  /// \short Helper function to set Euclid options using a command line
419  /// like array.
420  //=====================================================================
421  void euclid_settings_helper(const bool &use_block_jacobi,
422  const bool &use_row_scaling,
423  const bool &use_ilut,
424  const int &level,
425  const double &drop_tol,
426  const int &print_level,
427  HYPRE_Solver &euclid_object)
428  {
429  // Easier to use C-arrays rather than std::strings because Euclid takes
430  // char** as argument and because C++ doesn't provide decent number to
431  // std::string conversion functions.
432 
433  int n_args = 0;
434  const char* args[22];
435 
436  // first argument is empty string
437  args[n_args++] = "";
438 
439  // switch on/off Block Jacobi ILU
440  args[n_args++] = "-bj";
441  if (use_block_jacobi)
442  {
443  args[n_args++] = "1";
444  }
445  else
446  {
447  args[n_args++] = "0";
448  }
449 
450  // switch on/off row scaling
451  args[n_args++] = "-rowScale";
452  if (use_row_scaling)
453  {
454  args[n_args++] = "1";
455  }
456  else
457  {
458  args[n_args++] = "0";
459  }
460 
461  // set level for ILU(k)
462  args[n_args++] = "-level";
463  char level_value[10];
464  sprintf(level_value,"%d",level);
465  args[n_args++] = level_value;
466 
467  // // set drop tol for ILU(k) factorization
468  // args[n_args++] = "-sparseA";
469  // char droptol[20];
470  // sprintf(droptol,"%f",drop_tol);
471  // args[n_args++] = droptol;
472 
473  // // set ILUT factorization if required
474  // if (use_ilut)
475  // {
476  // args[n_args++] = "-ilut";
477  // args[n_args++] = droptol;
478  // }
479 
480  // set printing of Euclid data
481  if (print_level == 0)
482  {
483  args[n_args++] = "-eu_stats";
484  args[n_args++] = "0";
485  args[n_args++] = "-eu_mem";
486  args[n_args++] = "0";
487  }
488  if (print_level == 1)
489  {
490  args[n_args++] = "-eu_stats";
491  args[n_args++] = "1";
492  args[n_args++] = "-eu_mem";
493  args[n_args++] = "0";
494  }
495  if (print_level == 2)
496  {
497  args[n_args++] = "-eu_stats";
498  args[n_args++] = "1";
499  args[n_args++] = "-eu_mem";
500  args[n_args++] = "1";
501  }
502 
503  // set next entry in array to null
504  args[n_args] = 0;
505 
506  // Send the parameters
507  HYPRE_EuclidSetParams(euclid_object, n_args, const_cast<char**>(args));
508  }
509 
510  }
511 
512 
513 
514 ///////////////////////////////////////////////////////////////////////
515 ///////////////////////////////////////////////////////////////////////
516 // functions for HypreInterface class
517 ///////////////////////////////////////////////////////////////////////
518 ///////////////////////////////////////////////////////////////////////
519 
520 
521 
522 //=============================================================================
523 /// Helper function which creates a Hypre matrix from a CRDoubleMatrix
524 /// If OOMPH-LIB has been set up for MPI use, the Hypre matrix is
525 /// distributed over the available processors.
526 //=============================================================================
528  {
529 
530  // reset Hypre's global error flag
532 
533  // issue warning if the matrix is small compared to the number of processors
534  if ( unsigned(2*matrix_pt->distribution_pt()->communicator_pt()->nproc()) >
535  matrix_pt->nrow() )
536  {
537  oomph_info
538  << "Warning: HYPRE based solvers may fail if 2*number of processors "
539  << "is greater than the number of unknowns!" << std::endl;
540  }
541 
542  // store the distribution
543  // generate the Hypre matrix
545  Matrix_ij,
546  Matrix_par,
548 
549  // Output error messages if required
551  {
552  std::ostringstream message;
553  int err = HypreHelpers::check_HYPRE_error_flag(message);
554  if (err)
555  {
556  OomphLibWarning(message.str(),
557  "HypreSolver::hypre_matrix_setup()",
558  OOMPH_EXCEPTION_LOCATION);
559  }
560  }
561 
562  // delete CRDoubleMatrix if required
563  if (Delete_input_data)
564  {
565  matrix_pt->clear();
566  }
567  }
568 
569 
570 //=============================================================================
571 /// Sets up the solver data required for use in an oomph-lib
572 /// LinearSolver or Preconditioner, once the Hypre matrix has been
573 /// generated using hypre_matrix_setup(...).
574 //=============================================================================
576  {
577  // Store time
578  double t_start = TimingHelpers::timer();
579  double t_end = 0;
580 
581 
582  // reset Hypre's global error flag
584 
585  // create dummy Hypre vectors which are required for setup
586  HYPRE_IJVector dummy_sol_ij;
587  HYPRE_ParVector dummy_sol_par;
588  HYPRE_IJVector dummy_rhs_ij;
589  HYPRE_ParVector dummy_rhs_par;
591  dummy_sol_ij,
592  dummy_sol_par);
594  dummy_rhs_ij,
595  dummy_rhs_par);
596 
597  // Set up internal preconditioner for CG, GMRES or BiCGSTAB
598  // --------------------------------------------------------
599  if ( (Hypre_method>=CG) && (Hypre_method<=BiCGStab) )
600  {
601  // AMG preconditioner
603  {
604  // set up BoomerAMG
605  HYPRE_BoomerAMGCreate(&Preconditioner);
606  HYPRE_BoomerAMGSetPrintLevel(Preconditioner, AMG_print_level);
607  HYPRE_BoomerAMGSetMaxLevels(Preconditioner, AMG_max_levels);
608  HYPRE_BoomerAMGSetMaxIter(Preconditioner, 1);
609  HYPRE_BoomerAMGSetTol(Preconditioner, 0.0);
610  HYPRE_BoomerAMGSetCoarsenType(Preconditioner, AMG_coarsening);
611  HYPRE_BoomerAMGSetStrongThreshold(Preconditioner, AMG_strength);
612  HYPRE_BoomerAMGSetMaxRowSum(Preconditioner, AMG_max_row_sum);
613  HYPRE_BoomerAMGSetTruncFactor(Preconditioner, AMG_truncation);
614 
616  {
617  HYPRE_BoomerAMGSetRelaxType(Preconditioner, AMG_simple_smoother);
618  HYPRE_BoomerAMGSetNumSweeps(Preconditioner, AMG_smoother_iterations);
619 
620  // This one gives a memory leak
621  //double * relaxweight = new double[AMG_max_levels];
622 
623  // This is how they do it in a hypre demo code
624  double* relaxweight = hypre_CTAlloc(double, AMG_max_levels);
625 
626  for (unsigned i=0; i<AMG_max_levels; i++)
627  {
628  relaxweight[i] = AMG_damping;
629  }
630  HYPRE_BoomerAMGSetRelaxWeight(Preconditioner, relaxweight);
631  }
632  else
633  {
634 
635  HYPRE_BoomerAMGSetSmoothType(Preconditioner, AMG_complex_smoother);
636  HYPRE_BoomerAMGSetSmoothNumLevels(Preconditioner, AMG_max_levels);
637  HYPRE_BoomerAMGSetSmoothNumSweeps(Preconditioner,
639 
640  // If we are using Euclid then set up additional Euclid only options
641  if(AMG_complex_smoother == 9)
642  {
650  }
651  }
652 
654  }
655 
656  // Euclid preconditioner
658  {
659 #ifdef OOMPH_HAS_MPI
660  HYPRE_EuclidCreate(Hypre_distribution_pt->communicator_pt()->mpi_comm(),
661  &Preconditioner);
662 #else
663  HYPRE_EuclidCreate(MPI_COMM_WORLD, &Preconditioner);
664 #endif
665 
666  // Set parameters
670 
672  }
673 
674  // ParaSails preconditioner
676  {
677 #ifdef OOMPH_HAS_MPI
678  HYPRE_ParaSailsCreate
680 #else
681  HYPRE_ParaSailsCreate(MPI_COMM_WORLD, &Preconditioner);
682 #endif
683  HYPRE_ParaSailsSetSym(Preconditioner, ParaSails_symmetry);
684  HYPRE_ParaSailsSetParams(Preconditioner,
687  HYPRE_ParaSailsSetFilter(Preconditioner, ParaSails_filter);
689  }
690 
691  // check error flag
693  {
694  std::ostringstream message;
695  int err = HypreHelpers::check_HYPRE_error_flag(message);
696  if (err)
697  {
698  OomphLibWarning(message.str(),
699  "HypreSolver::hypre_setup()",
700  OOMPH_EXCEPTION_LOCATION);
701  }
702  }
703  } // end of setting up internal preconditioner
704 
705 
706  // set up solver
707  // -------------
708  t_start = TimingHelpers::timer();
709 
710  // AMG solver
711  if (Hypre_method==BoomerAMG)
712  {
713  if (Output_info)
714  {
715  oomph_info << "Setting up BoomerAMG, ";
716  }
717 
718  // set up BoomerAMG
719  HYPRE_BoomerAMGCreate(&Solver);
720  HYPRE_BoomerAMGSetPrintLevel(Solver, AMG_print_level);
721  HYPRE_BoomerAMGSetMaxLevels(Solver, AMG_max_levels);
722  HYPRE_BoomerAMGSetMaxIter(Solver, Max_iter);
723  HYPRE_BoomerAMGSetTol(Solver, Tolerance);
724  HYPRE_BoomerAMGSetCoarsenType(Solver, AMG_coarsening);
725  HYPRE_BoomerAMGSetStrongThreshold(Solver, AMG_strength);
726  HYPRE_BoomerAMGSetMaxRowSum(Solver, AMG_max_row_sum);
727  HYPRE_BoomerAMGSetTruncFactor(Solver, AMG_truncation);
728 
730  {
731  HYPRE_BoomerAMGSetRelaxType(Solver, AMG_simple_smoother);
732  HYPRE_BoomerAMGSetNumSweeps(Solver, AMG_smoother_iterations);
733 
734  // This one gives a memory leak
735  //double * relaxweight = new double[AMG_max_levels];
736 
737  // This is how they do it in a hypre demo code
738  double* relaxweight = hypre_CTAlloc(double, AMG_max_levels);
739 
740  for (unsigned i=0; i<AMG_max_levels; i++)
741  {
742  relaxweight[i] = AMG_damping;
743  }
744  HYPRE_BoomerAMGSetRelaxWeight(Solver, relaxweight);
745  }
746  else
747  {
748  HYPRE_BoomerAMGSetSmoothType(Solver, AMG_complex_smoother);
749  HYPRE_BoomerAMGSetSmoothNumLevels(Solver, AMG_max_levels);
750  HYPRE_BoomerAMGSetSmoothNumSweeps(Solver, AMG_smoother_iterations);
751 
752  /* Other settings
753  * 6 & Schwarz smoothers & HYPRE_BoomerAMGSetDomainType, HYPRE_BoomerAMGSetOverlap, \\
754  * & & HYPRE_BoomerAMGSetVariant, HYPRE_BoomerAMGSetSchwarzRlxWeight \\
755  * 7 & Pilut & HYPRE_BoomerAMGSetDropTol, HYPRE_BoomerAMGSetMaxNzPerRow \\
756  * 8 & ParaSails & HYPRE_BoomerAMGSetSym, HYPRE_BoomerAMGSetLevel, \\
757  * & & HYPRE_BoomerAMGSetFilter, HYPRE_BoomerAMGSetThreshold \\
758  * 9 & Euclid & HYPRE_BoomerAMGSetEuclidFile \\
759  */
760 
761  // If we are using Euclid then set up additional Euclid only options
762  if(AMG_complex_smoother == 9)
763  {
771  }
772 
773  // Add any others here as required...
774  }
775 
776 // MemoryUsage::doc_memory_usage("before amg setup [solver]");
777 // MemoryUsage::insert_comment_to_continous_top("BEFORE AMG SETUP [SOLVER]");
778 
779  HYPRE_BoomerAMGSetup(Solver,
780  Matrix_par,
781  dummy_rhs_par,
782  dummy_sol_par);
783 
784 // MemoryUsage::doc_memory_usage("after amg setup [solver]");
785 // MemoryUsage::insert_comment_to_continous_top("AFTER AMG SETUP [SOLVER]");
786 
788  }
789 
790  // Euclid solver
791  else if (Hypre_method==Euclid)
792  {
793  if (Output_info)
794  {
795  oomph_info << "Setting up Euclid, ";
796  }
797 #ifdef OOMPH_HAS_MPI
798  HYPRE_EuclidCreate(Hypre_distribution_pt->communicator_pt()->mpi_comm(),
799  &Solver);
800 #else
801  HYPRE_EuclidCreate(MPI_COMM_WORLD, &Solver);
802 #endif
803 
804  // Set parameters
808 
809  HYPRE_EuclidSetup(Solver,
810  Matrix_par,
811  dummy_rhs_par,
812  dummy_sol_par);
814  }
815 
816  // ParaSails preconditioner
817  else if (Hypre_method==ParaSails)
818  {
819  if (Output_info)
820  {
821  oomph_info << "Setting up ParaSails, ";
822  }
823 #ifdef OOMPH_HAS_MPI
824  HYPRE_ParaSailsCreate(Hypre_distribution_pt->communicator_pt()->mpi_comm(),
825  &Solver);
826 #else
827  HYPRE_ParaSailsCreate(MPI_COMM_WORLD, &Solver);
828 #endif
829  HYPRE_ParaSailsSetSym(Solver, ParaSails_symmetry);
830  HYPRE_ParaSailsSetParams(Solver,
833  HYPRE_ParaSailsSetFilter(Solver, ParaSails_filter);
834 
835  HYPRE_ParaSailsSetup(Solver,
836  Matrix_par,
837  dummy_rhs_par,
838  dummy_sol_par);
840  }
841 
842  // CG solver
843  else if (Hypre_method==CG)
844  {
845  if (Output_info)
846  {
847  oomph_info << "Setting up CG, ";
848  }
849 
850 #ifdef OOMPH_HAS_MPI
851  HYPRE_ParCSRPCGCreate(Hypre_distribution_pt->communicator_pt()->mpi_comm(),
852  &Solver);
853 #else
854  HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &Solver);
855 #endif
856  HYPRE_PCGSetTol(Solver, Tolerance);
857  HYPRE_PCGSetLogging(Solver, 0);
858  HYPRE_PCGSetPrintLevel(Solver, Krylov_print_level);
859  HYPRE_PCGSetMaxIter(Solver, Max_iter);
860 
861  // set preconditioner
863  {
864  if (Output_info)
865  {
866  oomph_info << " with BoomerAMG preconditioner, ";
867  }
868 
869  HYPRE_PCGSetPrecond(Solver,
870  (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSolve,
871  (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSetup,
873  }
874  else if (Internal_preconditioner==Euclid) // Euclid
875  {
876  if (Output_info)
877  {
878  oomph_info << " with Euclid ILU preconditioner, ";
879  }
880 
881  HYPRE_PCGSetPrecond(Solver,
882  (HYPRE_PtrToSolverFcn) HYPRE_EuclidSolve,
883  (HYPRE_PtrToSolverFcn) HYPRE_EuclidSetup,
885  }
886  else if (Internal_preconditioner==ParaSails) // ParaSails
887  {
888  if (Output_info)
889  {
890  oomph_info << " with ParaSails approximate inverse preconditioner, ";
891  }
892 
893  HYPRE_PCGSetPrecond(Solver,
894  (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSolve,
895  (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSetup,
897  }
898  else
899  {
900  if (Output_info)
901  {
902  oomph_info << " with no preconditioner";
903  }
904  }
905 
906 
907  HYPRE_PCGSetup(Solver,
908  (HYPRE_Matrix) Matrix_par,
909  (HYPRE_Vector) dummy_rhs_par,
910  (HYPRE_Vector) dummy_sol_par);
911 
912 
914  }
915 
916  // GMRES solver
917  else if (Hypre_method==GMRES)
918  {
919  if (Output_info)
920  {
921  oomph_info << "Setting up GMRES";
922  }
923 
924 #ifdef OOMPH_HAS_MPI
925  HYPRE_ParCSRGMRESCreate
926  (Hypre_distribution_pt->communicator_pt()->mpi_comm(),&Solver);
927 #else
928  HYPRE_ParCSRGMRESCreate(MPI_COMM_WORLD,&Solver);
929 #endif
930  HYPRE_GMRESSetTol(Solver, Tolerance);
931  HYPRE_GMRESSetKDim(Solver, Max_iter);
932  HYPRE_GMRESSetLogging(Solver, 0);
933  HYPRE_GMRESSetPrintLevel(Solver, Krylov_print_level);
934  HYPRE_GMRESSetMaxIter(Solver, Max_iter);
935 
936  // set preconditioner
938  {
939  if (Output_info)
940  {
941  oomph_info << " with BoomerAMG preconditioner, ";
942  }
943 
944  HYPRE_GMRESSetPrecond(Solver,
945  (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSolve,
946  (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSetup,
948  }
949  else if (Internal_preconditioner==Euclid) // Euclid
950  {
951  if (Output_info)
952  {
953  oomph_info << " with Euclid ILU preconditioner, ";
954  }
955 
956  HYPRE_GMRESSetPrecond(Solver,
957  (HYPRE_PtrToSolverFcn) HYPRE_EuclidSolve,
958  (HYPRE_PtrToSolverFcn) HYPRE_EuclidSetup,
960  }
961  else if (Internal_preconditioner==ParaSails) // ParaSails
962  {
963  if (Output_info)
964  {
965  oomph_info << " with ParaSails approximate inverse preconditioner, ";
966  }
967 
968  HYPRE_GMRESSetPrecond(Solver,
969  (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSolve,
970  (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSetup,
972  }
973  else
974  {
975  if (Output_info)
976  {
977  oomph_info << " with no preconditioner";
978  }
979  }
980 
981  HYPRE_GMRESSetup(Solver,
982  (HYPRE_Matrix) Matrix_par,
983  (HYPRE_Vector) dummy_rhs_par,
984  (HYPRE_Vector) dummy_sol_par);
985 
987  }
988 
989  // BiCGStab solver
990  else if (Hypre_method==BiCGStab)
991  {
992  if (Output_info)
993  {
994  oomph_info << "Setting up BiCGStab";
995  }
996 #ifdef OOMPH_HAS_MPI
997  HYPRE_ParCSRBiCGSTABCreate
998  (Hypre_distribution_pt->communicator_pt()->mpi_comm(), &Solver);
999 #else
1000  HYPRE_ParCSRBiCGSTABCreate(MPI_COMM_WORLD, &Solver);
1001 #endif
1002  HYPRE_BiCGSTABSetTol(Solver, Tolerance);
1003  HYPRE_BiCGSTABSetLogging(Solver, 0);
1004  HYPRE_BiCGSTABSetPrintLevel(Solver, Krylov_print_level);
1005  HYPRE_BiCGSTABSetMaxIter(Solver, Max_iter);
1006 
1007  // set preconditioner
1008  if (Internal_preconditioner==BoomerAMG) // AMG
1009  {
1010  if (Output_info)
1011  {
1012  oomph_info << " with BoomerAMG preconditioner, ";
1013  }
1014 
1015  HYPRE_BiCGSTABSetPrecond(Solver,
1016  (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSolve,
1017  (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSetup,
1018  Preconditioner);
1019  }
1020  else if (Internal_preconditioner==Euclid) // Euclid
1021  {
1022  if (Output_info)
1023  {
1024  oomph_info << " with Euclid ILU preconditioner, ";
1025  }
1026 
1027  HYPRE_BiCGSTABSetPrecond(Solver,
1028  (HYPRE_PtrToSolverFcn) HYPRE_EuclidSolve,
1029  (HYPRE_PtrToSolverFcn) HYPRE_EuclidSetup,
1030  Preconditioner);
1031  }
1032  else if (Internal_preconditioner==ParaSails) // ParaSails
1033  {
1034  if (Output_info)
1035  {
1036  oomph_info << " with ParaSails approximate inverse preconditioner, ";
1037  }
1038 
1039  HYPRE_BiCGSTABSetPrecond(Solver,
1040  (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSolve,
1041  (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSetup,
1042  Preconditioner);
1043  }
1044  else
1045  {
1046  if (Output_info)
1047  {
1048  oomph_info << " with no preconditioner, ";
1049  }
1050  }
1051 
1052  HYPRE_BiCGSTABSetup(Solver,
1053  (HYPRE_Matrix) Matrix_par,
1054  (HYPRE_Vector) dummy_rhs_par,
1055  (HYPRE_Vector) dummy_sol_par);
1056 
1058  }
1059 
1060  // no solver exists for this value of Solver flag
1061  else
1062  {
1063  std::ostringstream error_message;
1064  error_message << "Solver has been set to an invalid value. "
1065  << "current value=" << Solver;
1066  throw OomphLibError(error_message.str(),
1067  OOMPH_CURRENT_FUNCTION,
1068  OOMPH_EXCEPTION_LOCATION);
1069 
1070  }
1071 
1072  t_end = TimingHelpers::timer();
1073  double solver_setup_time = t_end-t_start;
1074 
1075  // destroy dummy hypre vectors
1076  HYPRE_IJVectorDestroy(dummy_sol_ij);
1077  HYPRE_IJVectorDestroy(dummy_rhs_ij);
1078 
1079  // check error flag
1081  {
1082  std::ostringstream message;
1083  int err = HypreHelpers::check_HYPRE_error_flag(message);
1084  if (err)
1085  {
1086  OomphLibWarning(message.str(),
1087  "HypreSolver::hypre_solver_setup()",
1088  OOMPH_EXCEPTION_LOCATION);
1089  }
1090  }
1091 
1092  // output times
1093  if (Output_info)
1094  {
1095  oomph_info << "time for setup [s] : "
1096  << solver_setup_time << std::endl;
1097  }
1098  }
1099 
1100 
1101 
1102 //===================================================================
1103 /// Helper function performs a solve if solver data has been set
1104 /// up using hypre_solver_setup(...).
1105 //====================================================================
1107  DoubleVector &solution)
1108  {
1109  // Record time
1110  double t_start = TimingHelpers::timer();
1111 
1112  // Set up hypre vectors
1113  // --------------------
1114 
1115  // Hypre vector for rhs values
1116  HYPRE_IJVector rhs_ij;
1117  HYPRE_ParVector rhs_par;
1118 
1119  // Hypre vector for solution
1120  HYPRE_IJVector solution_ij;
1121  HYPRE_ParVector solution_par;
1122 
1123  // set up rhs_values and vec_indices
1126  rhs_ij,
1127  rhs_par);
1128 
1130  solution_ij,
1131  solution_par);
1132 
1133  // check error flag
1135  {
1136  std::ostringstream message;
1137  int err = HypreHelpers::check_HYPRE_error_flag(message);
1138  if (err)
1139  {
1140  OomphLibWarning(message.str(),
1141  "HypreSolver::hypre_solve()",
1142  OOMPH_EXCEPTION_LOCATION);
1143  }
1144  }
1145 
1146  // solve
1147  // -----
1148 
1149  // for solver stats
1150  int iterations=0;
1151  double norm=0;
1152 
1153  // Get the norm of rhs
1154  const double rhs_norm = rhs.norm();
1155  bool do_solving = false;
1156  if (rhs_norm > 0.0)
1157  {
1158  do_solving = true;
1159  }
1160 
1161 #ifdef OOMPH_HAS_MPI
1162  // We need to check whether any processor requires to solve, if that
1163  // is the case then do the solving
1165  {
1166  if (MPI_Helpers::communicator_pt()->nproc()>1)
1167  {
1168  unsigned this_processor_do_solving = 0;
1169  unsigned all_processors_do_solving = 0;
1170  if (do_solving)
1171  {
1172  this_processor_do_solving = 1;
1173  }
1174  // Get the communicator
1176  // Communicate with all procesoors
1177  MPI_Allreduce(&this_processor_do_solving,&all_processors_do_solving,
1178  1, MPI_UNSIGNED,MPI_SUM,comm_pt->mpi_comm());
1179  if (all_processors_do_solving > 0)
1180  {
1181  do_solving = true;
1182  }
1183  }
1184  }
1185 #endif
1186 
1187  if (do_solving)
1188  {
1190  {
1191  HYPRE_BoomerAMGSolve(Solver, Matrix_par, rhs_par, solution_par);
1192  HYPRE_BoomerAMGGetNumIterations(Solver, &iterations);
1193  HYPRE_BoomerAMGGetFinalRelativeResidualNorm(Solver, &norm);
1194  }
1195  else if (Existing_solver==CG)
1196  {
1197  HYPRE_PCGSolve(Solver,
1198  (HYPRE_Matrix) Matrix_par,
1199  (HYPRE_Vector) rhs_par,
1200  (HYPRE_Vector) solution_par);
1201  HYPRE_PCGGetNumIterations(Solver, &iterations);
1202  HYPRE_PCGGetFinalRelativeResidualNorm(Solver, &norm);
1203  }
1204  else if (Existing_solver==GMRES)
1205  {
1206  HYPRE_GMRESSolve(Solver,
1207  (HYPRE_Matrix) Matrix_par,
1208  (HYPRE_Vector) rhs_par,
1209  (HYPRE_Vector) solution_par);
1210  HYPRE_GMRESGetNumIterations(Solver, &iterations);
1211  HYPRE_GMRESGetFinalRelativeResidualNorm(Solver, &norm);
1212  }
1213  else if (Existing_solver==BiCGStab)
1214  {
1215  HYPRE_BiCGSTABSolve(Solver,
1216  (HYPRE_Matrix) Matrix_par,
1217  (HYPRE_Vector) rhs_par,
1218  (HYPRE_Vector) solution_par);
1219  HYPRE_BiCGSTABGetNumIterations(Solver, &iterations);
1220  HYPRE_BiCGSTABGetFinalRelativeResidualNorm(Solver, &norm);
1221  }
1222  else if (Existing_solver==Euclid)
1223  {
1224  HYPRE_EuclidSolve(Solver, Matrix_par, rhs_par, solution_par);
1225  }
1226  else if (Existing_solver==ParaSails)
1227  {
1228  HYPRE_ParaSailsSolve(Solver, Matrix_par, rhs_par, solution_par);
1229  }
1230 
1231  // output any error message
1233  {
1234  std::ostringstream message;
1235  int err = HypreHelpers::check_HYPRE_error_flag(message);
1236  if (err)
1237  {
1238  OomphLibWarning(message.str(),
1239  "HypreSolver::hypre_solve()",
1240  OOMPH_EXCEPTION_LOCATION);
1241  }
1242  }
1243 
1244  } // if (do_solving)
1245 
1246  // Copy result to solution
1247  unsigned nrow_local = Hypre_distribution_pt->nrow_local();
1248  unsigned first_row = Hypre_distribution_pt->first_row();
1249  int* indices = new int[nrow_local];
1250  for (unsigned i = 0; i < nrow_local; i++)
1251  {
1252  indices[i] = first_row + i;
1253  }
1254  LinearAlgebraDistribution* soln_dist_pt;
1255  if (solution.built())
1256  {
1257  soln_dist_pt = new LinearAlgebraDistribution(solution.distribution_pt());
1258  }
1259  else
1260  {
1261  soln_dist_pt = new LinearAlgebraDistribution(rhs.distribution_pt());
1262  }
1263  solution.build(Hypre_distribution_pt,0.0);
1264  HYPRE_IJVectorGetValues(solution_ij,
1265  nrow_local,
1266  indices,
1267  solution.values_pt());
1268  solution.redistribute(soln_dist_pt);
1269  delete[] indices;
1270  delete soln_dist_pt;
1271 
1272  // output any error message
1274  {
1275  std::ostringstream message;
1276  int err = HypreHelpers::check_HYPRE_error_flag(message);
1277  if (err)
1278  {
1279  OomphLibWarning(message.str(),
1280  "HypreSolver::hypre_solve()",
1281  OOMPH_EXCEPTION_LOCATION);
1282  }
1283  }
1284 
1285  // deallocation
1286  HYPRE_IJVectorDestroy(solution_ij);
1287  HYPRE_IJVectorDestroy(rhs_ij);
1288 
1289  // Record time
1290  double solve_time=0;
1291  if (Output_info)
1292  {
1293  double t_end = TimingHelpers::timer();
1294  solve_time = t_end-t_start;
1295  }
1296 
1297  // output timings and info
1298  if (Output_info)
1299  {
1300  oomph_info << "Time for HYPRE solve [s] : "
1301  << solve_time << std::endl;
1302  }
1303 
1304  // for iterative solvers output iterations and final norm
1305  if ((Hypre_method>=CG) && (Hypre_method<=BoomerAMG))
1306  {
1307  if (iterations>1)
1308  {
1309  if (Output_info) oomph_info << "Number of iterations : "
1310  << iterations << std::endl;
1311  if (Output_info) oomph_info << "Final Relative Residual Norm : "
1312  << norm << std::endl;
1313  }
1314  }
1315  }
1316 
1317 
1318 //===================================================================
1319 /// hypre_clean_up_memory() deletes any existing Hypre solver and
1320 /// Hypre matrix
1321 //====================================================================
1323  {
1324  // is there an existing solver
1325  if (Existing_solver != None)
1326  {
1327 
1328  // delete matrix
1329  HYPRE_IJMatrixDestroy(Matrix_ij);
1330 
1331  // delete solver
1333  {
1334  HYPRE_BoomerAMGDestroy(Solver);
1335  }
1336  else if (Existing_solver==CG)
1337  {
1338  HYPRE_ParCSRPCGDestroy(Solver);
1339  }
1340  else if (Existing_solver==GMRES)
1341  {
1342  HYPRE_ParCSRGMRESDestroy(Solver);
1343  }
1344  else if (Existing_solver==BiCGStab)
1345  {
1346  HYPRE_ParCSRBiCGSTABDestroy(Solver);
1347  }
1348  else if (Existing_solver==Euclid)
1349  {
1350  HYPRE_EuclidDestroy(Solver);
1351  }
1352  else if (Existing_solver==ParaSails)
1353  {
1354  HYPRE_ParaSailsDestroy(Solver);
1355  }
1357 
1358  // delete preconditioner
1360  {
1361  HYPRE_BoomerAMGDestroy(Preconditioner);
1362  }
1363  else if (Existing_preconditioner==Euclid)
1364  {
1365  HYPRE_EuclidDestroy(Preconditioner);
1366  }
1368  {
1369  HYPRE_ParaSailsDestroy(Preconditioner);
1370  }
1372 
1373  // check error flag
1375  {
1376  std::ostringstream message;
1377  int err = HypreHelpers::check_HYPRE_error_flag(message);
1378  if (err)
1379  {
1380  OomphLibWarning(message.str(),
1381  "HypreSolver::clean_up_memory()",
1382  OOMPH_EXCEPTION_LOCATION);
1383  }
1384  }
1385  }
1386  }
1387 
1388 
1389 ///////////////////////////////////////////////////////////////////////
1390 ///////////////////////////////////////////////////////////////////////
1391 // functions for HypreSolver class
1392 ///////////////////////////////////////////////////////////////////////
1393 ///////////////////////////////////////////////////////////////////////
1394 
1395 
1396 //===================================================================
1397 /// Problem-based solve function to generate the Jacobian matrix and
1398 /// residual vector and use HypreInterface::hypre_solver_setup(...)
1399 /// and HypreInterface::hypre_solve(...) to solve the linear system.
1400 /// This function will delete any existing data.
1401 /// Note: the returned solution vector is NOT distributed, i.e. all
1402 /// processors hold all values because this is what the Newton solver
1403 /// requires.
1404 //====================================================================
1405  void HypreSolver::solve(Problem* const &problem_pt,
1406  DoubleVector &solution)
1407  {
1408  double t_start = TimingHelpers::timer();
1409 
1410  // Set Output_time flag for HypreInterface
1412 
1413  // Delete any existing solver data
1414  clean_up_memory();
1415 
1416  // Set flag to allow deletion of the oomphlib Jacobian matrix
1417  // (we're in control)
1418  Delete_input_data = true;
1419 
1420  //Get oomph-lib Jacobian matrix and residual vector
1421  DoubleVector residual;
1422  CRDoubleMatrix* matrix = new CRDoubleMatrix;
1423  problem_pt->get_jacobian(residual,*matrix);
1424 
1425  // Output times
1426  if(Doc_time)
1427  {
1428  oomph_info << "Time to generate Jacobian and residual [s] : ";
1429  double t_end = TimingHelpers::timer();
1430  oomph_info << t_end-t_start << std::endl;
1431  }
1432 
1433  // generate hypre matrix
1434  hypre_matrix_setup(matrix);
1435 
1436  // call hypre_solver_setup to generate linear solver data
1438 
1439  // perform hypre_solve
1440  hypre_solve(residual, solution);
1441 
1442  // delete solver data if required
1443  if (!Enable_resolve)
1444  {
1445  clean_up_memory();
1446  }
1447  }
1448 
1449 
1450 
1451 //===================================================================
1452 /// Uses HypreInterface::hypre_solve(...) to solve the linear system
1453 /// for a CRDoubleMatrix or a DistributedCRDoubleMatrix.
1454 /// In the latter case, the rhs needs to be distributed too,
1455 /// i.e. the length of the rhs vector must be equal to the number of
1456 /// rows stored locally.
1457 /// Note: the returned solution vector is never distributed, i.e. all
1458 /// processors hold all values
1459 //====================================================================
1460  void HypreSolver::solve(DoubleMatrixBase* const &matrix_pt,
1461  const DoubleVector &rhs,
1462  DoubleVector &solution)
1463  {
1464 #ifdef PARANOID
1465  // check the matrix is square
1466  if ( matrix_pt->nrow() != matrix_pt->ncol() )
1467  {
1468  std::ostringstream error_message;
1469  error_message << "HypreSolver require a square matrix. "
1470  << "Matrix is " << matrix_pt->nrow()
1471  << " by " << matrix_pt->ncol() << std::endl;
1472  throw OomphLibError(error_message.str(),
1473  OOMPH_CURRENT_FUNCTION,
1474  OOMPH_EXCEPTION_LOCATION);
1475  }
1476 #endif
1477 
1478  // Set Output_time flag for HypreInterface
1480 
1481  // Clean up existing solver data
1482  clean_up_memory();
1483 
1484  // Set flag to decide if oomphlib matrix can be deleted
1485  // (Recall that Delete_matrix defaults to false).
1487 
1488  // Try cast to a CRDoubleMatrix
1489  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
1490 
1491  // If cast successful set things up for a serial solve
1492  if (cr_matrix_pt)
1493  {
1494  // rebuild the distribution
1495  this->build_distribution(cr_matrix_pt->distribution_pt());
1496 
1497 #ifdef PARANOID
1498  // check that rhs has the same distribution as the matrix (now stored
1499  // in Distribution_pt)
1500  if (*this->distribution_pt() != *rhs.distribution_pt())
1501  {
1502  std::ostringstream error_message;
1503  error_message << "The distribution of the rhs vector and the matrix "
1504  << " must be the same" << std::endl;
1505  throw OomphLibError(error_message.str(),
1506  OOMPH_CURRENT_FUNCTION,
1507  OOMPH_EXCEPTION_LOCATION);
1508  }
1509  // if the solution is setup make sure it has the same distribution as
1510  // the matrix as well
1511  if (solution.built())
1512  {
1513  if (*this->distribution_pt() != *solution.distribution_pt())
1514  {
1515  std::ostringstream error_message;
1516  error_message << "The distribution of the solution vector is setup "
1517  << "there it must be the same as the matrix."
1518  << std::endl;
1519  throw OomphLibError(error_message.str(),
1520  OOMPH_CURRENT_FUNCTION,
1521  OOMPH_EXCEPTION_LOCATION);
1522  }
1523  }
1524 #endif
1525 
1526  hypre_matrix_setup(cr_matrix_pt);
1527  }
1528  else
1529  {
1530 #ifdef PARANOID
1531  std::ostringstream error_message;
1532  error_message << "HypreSolver only work with "
1533  << "CRDoubleMatrix matrices"
1534  << std::endl;
1535  throw OomphLibError(error_message.str(),
1536  OOMPH_CURRENT_FUNCTION,
1537  OOMPH_EXCEPTION_LOCATION);
1538 #endif
1539  }
1540 
1541  // call hypre_setup to generate Hypre matrix and linear solver data
1543 
1544  // perform hypre_solve
1545  hypre_solve(rhs, solution);
1546 
1547  // delete solver data if required
1548  if (!Enable_resolve)
1549  {
1550  clean_up_memory();
1551  }
1552  }
1553 
1554 
1555 //===================================================================
1556 /// Resolve performs a linear solve using current solver data (if
1557 /// such data exists).
1558 //====================================================================
1560  DoubleVector &solution)
1561  {
1562 #ifdef PARANOID
1563  // check solver data exists
1564  if (existing_solver()==None)
1565  {
1566  std::ostringstream error_message;
1567  error_message << "resolve(...) requires that solver data has been "
1568  << "set up by a previous call to solve(...) after "
1569  << "a call to enable_resolve()" << std::endl;
1570  throw OomphLibError(error_message.str(),
1571  OOMPH_CURRENT_FUNCTION,
1572  OOMPH_EXCEPTION_LOCATION);
1573  }
1574  // check that rhs has the same distribution as the matrix (now stored
1575  // in Distribution_pt)
1576  if (*this->distribution_pt() != *rhs.distribution_pt())
1577  {
1578  std::ostringstream error_message;
1579  error_message << "The distribution of the rhs vector and the matrix "
1580  << " must be the same" << std::endl;
1581  throw OomphLibError(error_message.str(),
1582  OOMPH_CURRENT_FUNCTION,
1583  OOMPH_EXCEPTION_LOCATION);
1584  }
1585  // if the solution is setup make sure it has the same distribution as
1586  // the matrix as well
1587  if (solution.built())
1588  {
1589  if (*this->distribution_pt() != *solution.distribution_pt())
1590  {
1591  std::ostringstream error_message;
1592  error_message << "The distribution of the solution vector is setup "
1593  << "there it must be the same as the matrix."
1594  << std::endl;
1595  throw OomphLibError(error_message.str(),
1596  OOMPH_CURRENT_FUNCTION,
1597  OOMPH_EXCEPTION_LOCATION);
1598  }
1599  }
1600 #endif
1601 
1602  // Set Output_info flag for HypreInterface
1604 
1605  // solve
1606  hypre_solve(rhs, solution);
1607 
1608  // Note: never delete solver data as the preconditioner is typically
1609  // called repeatedly.
1610  }
1611 
1612 
1613 //===================================================================
1614 /// clean_up_memory() deletes any existing solver data
1615 //====================================================================
1617  {
1619  }
1620 
1621 
1622 
1623 ///////////////////////////////////////////////////////////////////////
1624 ///////////////////////////////////////////////////////////////////////
1625 // functions for HyprePreconditioner class
1626 ///////////////////////////////////////////////////////////////////////
1627 ///////////////////////////////////////////////////////////////////////
1628 
1629 
1630 //=============================================================================
1631 /// \short Static double that accumulates the preconditioner
1632 /// solve time of all instantiations of this class. Reset
1633 /// it manually, e.g. after every Newton solve, using
1634  /// reset_cumulative_solve_times().
1635 //=============================================================================
1637 
1638 //=============================================================================
1639  /// \short map of static doubles that accumulates the preconditioner
1640  /// solve time of all instantiations of this class, labeled by
1641 /// context string. Reset
1642  /// it manually, e.g. after every Newton solve, using
1643  /// reset_cumulative_solve_times().
1644 //=============================================================================
1646 
1647 //=============================================================================
1648  /// \short Static unsigned that accumulates the number of preconditioner
1649  /// solves of all instantiations of this class. Reset
1650  /// it manually, e.g. after every Newton solve, using
1651  /// reset_cumulative_solve_times().
1652 //=============================================================================
1654 
1655 //=============================================================================
1656  /// \short Static unsigned that accumulates the number of preconditioner
1657  /// solves of all instantiations of this class, labeled by
1658  /// context string. Reset
1659  /// it manually, e.g. after every Newton solve, using
1660  /// reset_cumulative_solve_times().
1661 //=============================================================================
1662  std::map<std::string,unsigned>
1664 
1665 //=============================================================================
1666  /// \short Static unsigned that stores nrow for the most recent
1667  /// instantiations of this class, labeled by
1668  /// context string. Reset
1669  /// it manually, e.g. after every Newton solve, using
1670  /// reset_cumulative_solve_times().
1671 //=============================================================================
1672  std::map<std::string,unsigned> HyprePreconditioner::Context_based_nrow;
1673 
1674 //=============================================================================
1675  /// \short Report cumulative solve times of all instantiations of this
1676  /// class
1677 //=============================================================================
1679  {
1680  oomph_info << "\n\n=====================================================\n";
1681  oomph_info << "Cumulative HyprePreconditioner solve time "
1683  << " for " << Cumulative_npreconditioner_solve
1684  << " solves";
1685  if (Cumulative_npreconditioner_solve!=0)
1686  {
1687  oomph_info <<" ( "
1689  double(Cumulative_npreconditioner_solve) << " per solve )";
1690  }
1691  oomph_info << std::endl << std::endl;
1693  {
1694  oomph_info << "Breakdown by context: " << std::endl;
1695  for (std::map<std::string,double>::iterator it=
1698  {
1699  oomph_info << (*it).first << " " << (*it).second << " for "
1701  << " solves";
1703  {
1704  oomph_info << " ( "
1705  << (*it).second/
1707  << " per solve; "
1708  << (*it).second/
1710  double(Context_based_nrow[(*it).first])
1711  << " per solve per dof )";
1712  }
1713  oomph_info << std::endl;
1714  }
1715  }
1716  oomph_info << "\n=====================================================\n";
1717  oomph_info << std::endl;
1718  }
1719 
1720 //=============================================================================
1721  /// \short Reset cumulative solve times
1722 //=============================================================================
1724  {
1729  Context_based_nrow.clear();
1730  }
1731 
1732 
1733 //=============================================================================
1734 /// An interface to allow HypreSolver to be used as a Preconditioner
1735 /// for the oomph-lib IterativeLinearSolver class.
1736 /// Matrix has to be of type CRDoubleMatrix or DistributedCRDoubleMatrix.
1737 //=============================================================================
1739  {
1740  // Set Output_info flag for HypreInterface
1742 
1743 #ifdef PARANOID
1744  // check the matrix is square
1745  if ( matrix_pt()->nrow() != matrix_pt()->ncol() )
1746  {
1747  std::ostringstream error_message;
1748  error_message << "HyprePreconditioner require a square matrix. "
1749  << "Matrix is " << matrix_pt()->nrow()
1750  << " by " << matrix_pt()->ncol() << std::endl;
1751  throw OomphLibError(error_message.str(),
1752  OOMPH_CURRENT_FUNCTION,
1753  OOMPH_EXCEPTION_LOCATION);
1754  }
1755 #endif
1756 
1757  // clean up any previous solver data
1758  clean_up_memory();
1759 
1760  // set flag to decide if oomphlib matrix can be deleted
1761  // (Recall that Delete_matrix defaults to false).
1763 
1764  // Try casting to a CRDoubleMatrix
1765  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
1766 
1767  // If cast successful set things up for a serial solve
1768  if (cr_matrix_pt)
1769  {
1770  this->build_distribution(cr_matrix_pt->distribution_pt());
1771  hypre_matrix_setup(cr_matrix_pt);
1772  }
1773  else
1774  {
1775 #ifdef PARANOID
1776  std::ostringstream error_message;
1777  error_message << "HyprePreconditioner only work with "
1778  << "CRDoubleMatrix matrices"
1779  << std::endl;
1780  throw OomphLibError(error_message.str(),
1781  OOMPH_CURRENT_FUNCTION,
1782  OOMPH_EXCEPTION_LOCATION);
1783 #endif
1784  }
1785 
1786  if (Context_string!="")
1787  {
1788  oomph_info << "Setup of HyprePreconditioner in context \" "
1789  << Context_string << "\": nrow, nrow_local, nnz "
1790  << cr_matrix_pt->nrow() << " "
1791  << cr_matrix_pt->nrow_local() << " "
1792  << cr_matrix_pt->nnz() << std::endl;
1793  }
1794  Context_based_nrow[Context_string]=cr_matrix_pt->nrow();
1795 
1796  // call hypre_solver_setup
1798  }
1799 
1800 //===================================================================
1801 /// Preconditioner_solve uses a hypre solver to precondition vector r
1802 //====================================================================
1804  DoubleVector &z)
1805  {
1806 
1807  // Store time
1808  double t_start = TimingHelpers::timer();
1809 
1810 #ifdef PARANOID
1811  // check solver data exists
1812  if (existing_solver()==None)
1813  {
1814  std::ostringstream error_message;
1815  error_message << "preconditioner_solve(...) requires that data has "
1816  << "been set up using the function setup(...)" << std::endl;
1817  throw OomphLibError(error_message.str(),
1818  OOMPH_CURRENT_FUNCTION,
1819  OOMPH_EXCEPTION_LOCATION);
1820  }
1821  // check that rhs has the same distribution as the matrix (now stored
1822  // in Distribution_pt)
1823  if (*this->distribution_pt() != *r.distribution_pt())
1824  {
1825  std::ostringstream error_message;
1826  error_message << "The distribution of the rhs vector and the matrix "
1827  << " must be the same" << std::endl;
1828  throw OomphLibError(error_message.str(),
1829  OOMPH_CURRENT_FUNCTION,
1830  OOMPH_EXCEPTION_LOCATION);
1831  }
1832 
1833  // if the solution is setup make sure it has the same distribution as
1834  // the matrix as well
1835  if (z.built())
1836  {
1837  if (*this->distribution_pt() != *z.distribution_pt())
1838  {
1839  std::ostringstream error_message;
1840  error_message << "The distribution of the solution vector is setup "
1841  << "there it must be the same as the matrix."
1842  << std::endl;
1843  throw OomphLibError(error_message.str(),
1844  OOMPH_CURRENT_FUNCTION,
1845  OOMPH_EXCEPTION_LOCATION);
1846  }
1847  }
1848 #endif
1849 
1850  // Switch off any timings for the solve
1851  Output_info = false;
1852 
1853  // perform hypre_solve
1854  hypre_solve(r,z);
1855 
1856  // Add to cumulative solve time
1857  double t_end = TimingHelpers::timer();
1858  Cumulative_preconditioner_solve_time+=(t_end-t_start);
1860  My_cumulative_preconditioner_solve_time+=(t_end-t_start);
1861  if (Context_string!="")
1862  {
1865  }
1866  }
1867 
1868 
1869 
1870 //===================================================================
1871 /// clean_up_memory() deletes any existing Hypre solver and
1872 /// Hypre matrix
1873 //====================================================================
1875  {
1877  }
1878 
1879 }
bool Output_info
Flag is true to output info and results of timings.
Definition: hypre_solver.h:307
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1056
unsigned Krylov_print_level
Used to set the Hypre printing level for the Krylov subspace solvers.
Definition: hypre_solver.h:431
unsigned & amg_simple_smoother()
Access function to AMG_simple_smoother flag.
Definition: hypre_solver.h:888
int Euclid_level
Euclid level parameter for ILU(k) factorization.
Definition: hypre_solver.h:420
bool built() const
unsigned nrow_local() const
access function for the num of local rows on this processor.
unsigned Hypre_method
Hypre method flag. Valid values are specified in enumeration.
Definition: hypre_solver.h:316
bool Hypre_error_messages
Flag to determine if non-zero values of the Hypre error flag plus Hypre error messages are output to ...
Definition: hypre_solver.h:443
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero...
HYPRE_IJMatrix Matrix_ij
The Hypre_IJMatrix version of the matrix used in solve(...), resolve(...) or preconditioner_solve(...).
Definition: hypre_solver.h:472
HYPRE_Solver Preconditioner
The internal Hypre preconditioner used in conjunction with Solver. [This is a C structure!].
Definition: hypre_solver.h:484
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
bool AMGEuclidSmoother_use_block_jacobi
Definition: hypre_solver.h:280
void clean_up_memory()
Function deletes all solver data.
double & amg_damping()
Access function to AMG_damping parameter.
Definition: hypre_solver.h:916
double * values_pt()
access function to the underlying values
int hypre__global_error
cstr elem_len * i
Definition: cfortran.h:607
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
void euclid_settings_helper(const bool &use_block_jacobi, const bool &use_row_scaling, const bool &use_ilut, const int &level, const double &drop_tol, const int &print_level, HYPRE_Solver &euclid_object)
Helper function to set Euclid options using a command line like array.
bool AMG_using_simple_smoothing
Flag to determine whether simple smoothers (determined by the AMG_simple_smoother flag) or complex sm...
Definition: hypre_solver.h:340
unsigned nrow() const
access function to the number of global rows.
unsigned AMG_print_level
Used to set the Hypre printing level for AMG 0: no printout 1: print setup information 2: print solve...
Definition: hypre_solver.h:329
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned...
std::string Context_string
String can be used to provide context for annotation.
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:80
The Problem class.
Definition: problem.h:152
double AMGEuclidSmoother_drop_tol
Definition: hypre_solver.h:284
int ParaSails_nlevel
ParaSails nlevel parameter.
Definition: hypre_solver.h:399
void build(const OomphCommunicator *const comm_pt, const unsigned &first_row, const unsigned &nrow_local, const unsigned &nrow=0)
Sets the distribution. Takes first_row, nrow_local and nrow as arguments. If nrow is not provided or ...
double My_cumulative_preconditioner_solve_time
Private double that accumulates the preconditioner solve time of thi instantiation of this class...
double Euclid_droptol
Euclid drop tolerance for ILU(k) and ILUT factorization.
Definition: hypre_solver.h:408
double AMG_strength
Connection strength threshold parameter for BoomerAMG.
Definition: hypre_solver.h:371
unsigned & amg_coarsening()
Access function to AMG_coarsening flag.
Definition: hypre_solver.h:910
bool distributed() const
distribution is serial or distributed
static OomphCommunicator * communicator_pt()
access to the global oomph-lib communicator
OomphInfo oomph_info
void hypre_matrix_setup(CRDoubleMatrix *matrix_pt)
Function which sets values of First_global_row, Last_global_row and other partitioning data and creat...
static bool mpi_has_been_initialised()
return true if MPI has been initialised
The GMRES method.
unsigned Existing_preconditioner
Used to keep track of which preconditioner (if any) is currently stored.
Definition: hypre_solver.h:490
void set_amg_iterations(const unsigned &amg_iterations)
Function to set the number of times to apply BoomerAMG.
Definition: hypre_solver.h:872
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed ...
double * value()
Access to C-style value array.
Definition: matrices.h:1062
bool AMGEuclidSmoother_use_row_scaling
Definition: hypre_solver.h:281
unsigned Existing_solver
Used to keep track of which solver (if any) is currently stored.
Definition: hypre_solver.h:487
The conjugate gradient method.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned AMG_simple_smoother
Simple smoothing methods used in BoomerAMG. Relaxation types include: 0 = Jacobi 1 = Gauss-Seidel...
Definition: hypre_solver.h:353
unsigned Euclid_print_level
Flag to set the level of printing from Euclid when the Euclid destructor is called 0: no printing (de...
Definition: hypre_solver.h:427
double norm() const
compute the 2 norm of this vector
void solve(Problem *const &problem_pt, DoubleVector &solution)
Function which uses problem_pt's get_jacobian(...) function to generate a linear system which is then...
bool Euclid_rowScale
Flag to switch on Euclid row scaling.
Definition: hypre_solver.h:411
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
Definition: matrices.h:1068
double AMG_strength
Default for AMG strength (0.25 recommended for 2D problems; larger (0.5-0.75, say) for 3D...
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
void clear()
clear
Definition: matrices.cc:1677
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void amg_using_simple_smoothing()
Function to select use of 'simple' AMG smoothers as controlled by the flag AMG_simple_smoother.
Definition: hypre_solver.h:885
double AMG_truncation
Interpolation truncation factor for BoomerAMG.
Definition: hypre_solver.h:374
static unsigned Cumulative_npreconditioner_solve
Static unsigned that accumulates the number of preconditioner solves of all instantiations of this cl...
Definition: hypre_solver.h:823
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:84
void hypre_solver_setup()
Sets up the data required for to use as an oomph-lib LinearSolver or Preconditioner. This must be called after the Hypre matrix has been generated using hypre_matrix_setup(...).
unsigned AMGEuclidSmoother_print_level
Definition: hypre_solver.h:285
unsigned Max_iter
Maximum number of iterations used in solver.
Definition: hypre_solver.h:310
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:998
unsigned Internal_preconditioner
Preconditioner method flag used with Hypre's PCG, GMRES and BiCGStab in solve(...) or resolve(...
Definition: hypre_solver.h:322
static std::map< std::string, double > Context_based_cumulative_solve_time
Static double that accumulates the preconditioner solve time of all instantiations of this class...
Definition: hypre_solver.h:816
void create_HYPRE_Vector(const DoubleVector &oomph_vec, const LinearAlgebraDistribution *dist_pt, HYPRE_IJVector &hypre_ij_vector, HYPRE_ParVector &hypre_par_vector)
Helper function to create a HYPRE_IJVector and HYPRE_ParVector.
double AMG_damping
Damping factor for BoomerAMG smoothed Jacobi or hybrid SOR.
Definition: hypre_solver.h:368
void set_defaults_for_3D_poisson_problem(HyprePreconditioner *hypre_preconditioner_pt)
Set default parameters for use as preconditioner in 3D Poisson-type problem.
Definition: hypre_solver.cc:77
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Function applies solver to vector r for preconditioning. This requires a call to setup(...) first. Note: Hypre copies matrix data from oomph-lib's CRDoubleMatrix or DistributedCRDoubleMatrix into its own data structures, doubling the memory requirements for the matrix. As far as the Hypre solvers are concerned, the oomph-lib matrix is no longer required and could be deleted to save memory. However, this will not be the expected behaviour and therefore needs to be requested explicitly by the user with the delete_matrix() function.
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:992
unsigned AMG_smoother_iterations
The number of smoother iterations to apply.
Definition: hypre_solver.h:365
unsigned AMG_max_levels
Maximum number of levels used in AMG.
Definition: hypre_solver.h:332
unsigned AMG_coarsening
AMG coarsening strategy. Coarsening types include: 0 = CLJP (parallel coarsening using independent se...
Definition: hypre_solver.h:388
unsigned & hypre_method()
Access function to Hypre_method flag – specified via enumeration.
Definition: hypre_solver.h:862
unsigned existing_solver()
Function to return value of which solver (if any) is currently stored.
Definition: hypre_solver.h:274
void setup()
Function to set up a preconditioner for the linear system defined by matrix_pt. This function is requ...
double timer()
returns the time in seconds after some point in past
HYPRE_ParCSRMatrix Matrix_par
The Hypre_ParCSRMatrix version of the matrix used in solve(...), resolve(...) or preconditioner_solve...
Definition: hypre_solver.h:476
HYPRE_Solver Solver
The Hypre solver used in solve(...), resolve(...) or preconditioner_solve(...). [This is a C structur...
Definition: hypre_solver.h:480
bool Euclid_using_ILUT
Flag to determine if ILUT (if true) or ILU(k) is used in Euclid.
Definition: hypre_solver.h:414
unsigned AMG_complex_smoother
Complex smoothing methods used in BoomerAMG. Relaxation types are: 6 = Schwarz 7 = Pilut 8 = ParaSail...
Definition: hypre_solver.h:362
double ParaSails_thresh
ParaSails thresh parameter.
Definition: hypre_solver.h:402
bool built() const
access function to the Built flag - indicates whether the matrix has been build - i...
Definition: matrices.h:1177
unsigned AMGEuclidSmoother_level
Definition: hypre_solver.h:283
int ParaSails_symmetry
ParaSails symmetry flag, used to inform ParaSails of Symmetry of definitenss of problem and type of P...
Definition: hypre_solver.h:396
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
Definition: problem.cc:3851
The conjugate gradient method.
double & amg_strength()
Access function to AMG_strength.
Definition: hypre_solver.h:919
bool Delete_input_data
Internal flag which is true when hypre_setup or hypre_solve can delete input matrix.
Definition: hypre_solver.h:447
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
bool Euclid_using_BJ
Flag to determine if Block Jacobi is used instead of PILU.
Definition: hypre_solver.h:417
bool Delete_matrix
Hypre copies matrix data from oomph-lib's CRDoubleMatrix or DistributedCRDoubleMatrix into its own da...
static double Cumulative_preconditioner_solve_time
Static double that accumulates the preconditioner solve time of all instantiations of this class...
Definition: hypre_solver.h:809
void create_HYPRE_Matrix(CRDoubleMatrix *oomph_matrix, HYPRE_IJMatrix &hypre_ij_matrix, HYPRE_ParCSRMatrix &hypre_par_matrix, LinearAlgebraDistribution *dist_pt)
Helper function to create a serial HYPRE_IJMatrix and HYPRE_ParCSRMatrix from a CRDoubleMatrix.
void hypre_clean_up_memory()
Function deletes all solver data.
void hypre_solve(const DoubleVector &rhs, DoubleVector &solution)
Helper function performs a solve if any solver exists.
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
double AMG_truncation
AMG interpolation truncation factor.
static void reset_cumulative_solve_times()
Reset cumulative solve times.
void clean_up_memory()
Function deletes all solver data.
double AMG_max_row_sum
Parameter to identify diagonally dominant parts of the matrix in AMG.
Definition: hypre_solver.h:335
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
bool Delete_matrix
Hypre copies matrix data from oomph-lib's CRDoubleMatrix or DistributedCRDoubleMatrix into its own da...
Definition: hypre_solver.h:724
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
void set_defaults_for_navier_stokes_momentum_block(HyprePreconditioner *hypre_preconditioner_pt)
Set default parameters for use as preconditioner in for momentum block in Navier-Stokes problem...
Definition: hypre_solver.cc:90
static std::map< std::string, unsigned > Context_based_nrow
Static unsigned that stores nrow for the most recent instantiations of this class, labeled by context string. Reset it manually, e.g. after every Newton solve, using reset_cumulative_solve_times().
Definition: hypre_solver.h:838
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1050
void resolve(const DoubleVector &rhs, DoubleVector &solution)
Function to resolve a linear system using the existing solver data, allowing a solve with a new right...
unsigned AMG_coarsening
Default AMG coarsening strategy. Coarsening types include: 0 = CLJP (parallel coarsening using indepe...
LinearAlgebraDistribution * Hypre_distribution_pt
the distribution for this helpers-
Definition: hypre_solver.h:493
Abstract base class for matrices of doubles – adds abstract interfaces for solving, LU decomposition and multiplication by vectors.
Definition: matrices.h:275
double Tolerance
Tolerance used to terminate solver.
Definition: hypre_solver.h:313
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
static void report_cumulative_solve_times()
Report cumulative solve times of all instantiations of this class.
void set_defaults_for_2D_poisson_problem(HyprePreconditioner *hypre_preconditioner_pt)
Set default parameters for use as preconditioner in 2D Poisson-type problem.
Definition: hypre_solver.cc:50
static std::map< std::string, unsigned > Context_based_cumulative_npreconditioner_solve
Static unsigned that accumulates the number of preconditioner solves of all instantiations of this cl...
Definition: hypre_solver.h:831
double ParaSails_filter
ParaSails filter parameter.
Definition: hypre_solver.h:405
int check_HYPRE_error_flag(std::ostringstream &message)
Helper function to check the Hypre error flag, return the message associated with any error...
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57