oomph::HypreHelpers Namespace Reference

Helper functions for use with the Hypre library. More...

Functions
int	check_HYPRE_error_flag (std::ostringstream &message)
void	create_HYPRE_Vector (const DoubleVector &oomph_vec, const LinearAlgebraDistribution *dist_pt, HYPRE_IJVector &hypre_ij_vector, HYPRE_ParVector &hypre_par_vector)
void	create_HYPRE_Vector (const LinearAlgebraDistribution *dist_pt, HYPRE_IJVector &hypre_ij_vector, HYPRE_ParVector &hypre_par_vector)
void	create_HYPRE_Matrix (CRDoubleMatrix *oomph_matrix, HYPRE_IJMatrix &hypre_ij_matrix, HYPRE_ParCSRMatrix &hypre_par_matrix, LinearAlgebraDistribution *dist_pt)
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)

Variables
double	AMG_strength = 0.25
unsigned	AMG_coarsening = 6
double	AMG_truncation = 0.0
	AMG interpolation truncation factor. More...

Detailed Description

Helper functions for use with the Hypre library.

///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////

Function Documentation

check_HYPRE_error_flag()

int oomph::HypreHelpers::check_HYPRE_error_flag

(

std::ostringstream &

message

)

Helper function to check the Hypre error flag, return the message associated with any error, and reset the global error flag to zero. This function also returns the error value.

Helper function to check the Hypre error flag, return the message associated with any error, and reset the error flag to zero.

    {
      // get the Hypre error flag
      int err = HYPRE_GetError();
 
      // Tell us all about it...
      if (err)
      {
        oomph_info << "Hypre error flag=" << err << std::endl;
        char* error_message = new char[128];
        HYPRE_DescribeError(err, error_message);
        message << "WARNING: " << std::endl
                << "HYPRE error message: " << error_message << std::endl;
        delete[] error_message;
      }
 
      // reset Hypre's global error flag
      hypre__global_error = 0;
 
      return err;
    }

References hypre__global_error, and oomph::oomph_info.

Referenced by oomph::HypreInterface::hypre_clean_up_memory(), oomph::HypreInterface::hypre_matrix_setup(), oomph::HypreInterface::hypre_solve(), and oomph::HypreInterface::hypre_solver_setup().

create_HYPRE_Matrix()

void oomph::HypreHelpers::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 NOTE: dist_pt is rebuilt to match the distribution of the hypre solver which is not necassarily the same as the oomph lib matrix

Helper function to create a serial HYPRE_IJMatrix and HYPRE_ParCSRMatrix from a CRDoubleMatrix

    {
#ifdef PARANOID
      // check that the matrix is built
      if (!oomph_matrix->built())
      {
        std::ostringstream error_message;
        error_message << "The matrix has not been built";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      // check the matrix is square
      if (oomph_matrix->nrow() != oomph_matrix->ncol())
      {
        std::ostringstream error_message;
        error_message << "create_HYPRE_Matrix require a square matrix. "
                      << "Matrix is " << oomph_matrix->nrow() << " by "
                      << oomph_matrix->ncol() << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
#ifdef OOMPH_HAS_MPI
      // Trap the case when we have compiled with MPI,
      // but are running in serial
      if (!MPI_Helpers::mpi_has_been_initialised())
      {
        std::ostringstream error_stream;
        error_stream
          << "Oomph-lib has been compiled with MPI support and "
          << "you are using HYPRE.\n"
          << "For this combination of flags, MPI must be initialised.\n"
          << "Call MPI_Helpers::init() in the "
          << "main() function of your driver code\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // find number of rows/columns
      const unsigned nrow = int(oomph_matrix->nrow());
 
      // get pointers to the matrix
 
      // column indices of matrix
      const int* matrix_cols = oomph_matrix->column_index();
 
      // entries of matrix
      const double* matrix_vals = oomph_matrix->value();
 
      // row starts
      const int* matrix_row_start = oomph_matrix->row_start();
 
      // build the distribution
      if (oomph_matrix->distribution_pt()->distributed())
      {
        dist_pt->build(oomph_matrix->distribution_pt());
      }
      else
      {
        bool distributed = true;
        if (oomph_matrix->distribution_pt()->communicator_pt()->nproc() == 1)
        {
          distributed = false;
        }
        dist_pt->build(oomph_matrix->distribution_pt()->communicator_pt(),
                       nrow,
                       distributed);
      }
 
      // initialize hypre matrix
      unsigned lower = dist_pt->first_row();
      unsigned upper = lower + dist_pt->nrow_local() - 1;
 
#ifdef OOMPH_HAS_MPI
      HYPRE_IJMatrixCreate(dist_pt->communicator_pt()->mpi_comm(),
                           lower,
                           upper,
                           lower,
                           upper,
                           &hypre_ij_matrix);
#else
      HYPRE_IJMatrixCreate(
        MPI_COMM_WORLD, lower, upper, lower, upper, &hypre_ij_matrix);
#endif
      HYPRE_IJMatrixSetObjectType(hypre_ij_matrix, HYPRE_PARCSR);
      HYPRE_IJMatrixInitialize(hypre_ij_matrix);
 
      // set up a row map
      // and first row / nrow_local
      const unsigned hypre_nrow_local = dist_pt->nrow_local();
      const unsigned hypre_first_row = dist_pt->first_row();
      int* ncols_per_row = new int[hypre_nrow_local];
      int* row_map = new int[hypre_nrow_local];
      for (unsigned i = 0; i < hypre_nrow_local; i++)
      {
        unsigned j = i;
        if (!oomph_matrix->distributed() && dist_pt->distributed())
        {
          j += hypre_first_row;
        }
        ncols_per_row[i] = matrix_row_start[j + 1] - matrix_row_start[j];
        row_map[i] = hypre_first_row + i;
      }
 
      // put values in HYPRE matrix
      int local_start = 0;
      if (!oomph_matrix->distributed() && dist_pt->distributed())
      {
        local_start += matrix_row_start[hypre_first_row];
      }
 
 
      HYPRE_IJMatrixSetValues(hypre_ij_matrix,
                              hypre_nrow_local,
                              ncols_per_row,
                              row_map,
                              matrix_cols + local_start,
                              matrix_vals + local_start);
 
      // assemble matrix
      HYPRE_IJMatrixAssemble(hypre_ij_matrix); // hierher leak?
      HYPRE_IJMatrixGetObject(hypre_ij_matrix, (void**)&hypre_par_matrix);
 
      // tidy up memory
      delete[] ncols_per_row;
      delete[] row_map;
    }

References oomph::LinearAlgebraDistribution::build(), oomph::CRDoubleMatrix::built(), oomph::CRDoubleMatrix::column_index(), oomph::LinearAlgebraDistribution::communicator_pt(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearAlgebraDistribution::first_row(), i, int(), j, helpers::lower(), oomph::MPI_Helpers::mpi_has_been_initialised(), oomph::CRDoubleMatrix::ncol(), oomph::OomphCommunicator::nproc(), oomph::CRDoubleMatrix::nrow(), oomph::LinearAlgebraDistribution::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRDoubleMatrix::row_start(), and oomph::CRDoubleMatrix::value().

Referenced by oomph::HypreInterface::hypre_matrix_setup().

create_HYPRE_Vector() [1/2]

void oomph::HypreHelpers::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.

• If no MPI then serial vectors are created
• If MPI and serial input vector then distributed hypre vectors are created
• If MPI and distributed input vector the distributed output vectors are created.

    {
      // the lower and upper row of the vector on this processor
      unsigned lower = dist_pt->first_row();
      unsigned upper = dist_pt->first_row() + dist_pt->nrow_local() - 1;
 
      // number of local rows
      unsigned nrow_local = dist_pt->nrow_local();
 
      // initialize Hypre vector
#ifdef OOMPH_HAS_MPI
      HYPRE_IJVectorCreate(
        oomph_vec.distribution_pt()->communicator_pt()->mpi_comm(),
        lower,
        upper,
        &hypre_ij_vector);
#else
      HYPRE_IJVectorCreate(MPI_COMM_WORLD, lower, upper, &hypre_ij_vector);
#endif
      HYPRE_IJVectorSetObjectType(hypre_ij_vector, HYPRE_PARCSR);
      HYPRE_IJVectorInitialize(hypre_ij_vector);
 
      // set up array containing indices
      int* indices = new int[nrow_local];
      double* values = new double[nrow_local];
      const unsigned hypre_first_row = dist_pt->first_row();
      unsigned j = 0;
      if (!oomph_vec.distributed() && dist_pt->distributed())
      {
        j = hypre_first_row;
      }
      const double* o_pt = oomph_vec.values_pt();
      for (unsigned i = 0; i < nrow_local; i++)
      {
        indices[i] = hypre_first_row + i;
        values[i] = o_pt[j + i];
      }
 
      // insert values
      HYPRE_IJVectorSetValues(hypre_ij_vector, nrow_local, indices, values);
 
      // assemble vectors
      HYPRE_IJVectorAssemble(hypre_ij_vector);
      HYPRE_IJVectorGetObject(hypre_ij_vector, (void**)&hypre_par_vector);
 
      // clean up
      delete[] indices;
      delete[] values;
    }

References oomph::LinearAlgebraDistribution::communicator_pt(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearAlgebraDistribution::first_row(), i, j, helpers::lower(), oomph::LinearAlgebraDistribution::nrow_local(), and oomph::DoubleVector::values_pt().

Referenced by oomph::HypreInterface::hypre_solve(), and oomph::HypreInterface::hypre_solver_setup().

create_HYPRE_Vector() [2/2]

void oomph::HypreHelpers::create_HYPRE_Vector	(	const LinearAlgebraDistribution *	dist_pt,
		HYPRE_IJVector &	hypre_ij_vector,
		HYPRE_ParVector &	hypre_par_vector
	)

Helper function to create a HYPRE_IJVector and HYPRE_ParVector.

• If no MPI then serial vectors are created
• If MPI and serial input vector then distributed hypre vectors are created
• If MPI and distributed input vector the distributed output vectors are created.

Helper function to create an empty HYPRE_IJVector and HYPRE_ParVector.

• If no MPI then serial vectors are created
• If MPI and serial distribution then distributed hypre vectors are created
• If MPI and distributed input distribution the distributed output vectors are created.

    {
      // the lower and upper row of the vector on this processor
      unsigned lower = dist_pt->first_row();
      unsigned upper = dist_pt->first_row() + dist_pt->nrow_local() - 1;
 
      // initialize Hypre vector
#ifdef OOMPH_HAS_MPI
      HYPRE_IJVectorCreate(
        dist_pt->communicator_pt()->mpi_comm(), lower, upper, &hypre_ij_vector);
#else
      HYPRE_IJVectorCreate(MPI_COMM_WORLD, lower, upper, &hypre_ij_vector);
#endif
      HYPRE_IJVectorSetObjectType(hypre_ij_vector, HYPRE_PARCSR);
      HYPRE_IJVectorInitialize(hypre_ij_vector);
 
      // assemble vectors
      HYPRE_IJVectorAssemble(hypre_ij_vector);
      HYPRE_IJVectorGetObject(hypre_ij_vector, (void**)&hypre_par_vector);
    }

References oomph::LinearAlgebraDistribution::communicator_pt(), oomph::LinearAlgebraDistribution::first_row(), helpers::lower(), and oomph::LinearAlgebraDistribution::nrow_local().

euclid_settings_helper()

void oomph::HypreHelpers::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.

    {
      // Easier to use C-arrays rather than std::strings because Euclid takes
      // char** as argument and because C++ doesn't provide decent number to
      // std::string conversion functions.
 
      int n_args = 0;
      const char* args[22];
 
      // first argument is empty string
      args[n_args++] = "";
 
      // switch on/off Block Jacobi ILU
      args[n_args++] = "-bj";
      if (use_block_jacobi)
      {
        args[n_args++] = "1";
      }
      else
      {
        args[n_args++] = "0";
      }
 
      // switch on/off row scaling
      args[n_args++] = "-rowScale";
      if (use_row_scaling)
      {
        args[n_args++] = "1";
      }
      else
      {
        args[n_args++] = "0";
      }
 
      // set level for ILU(k)
      args[n_args++] = "-level";
      char level_value[10];
      sprintf(level_value, "%d", level);
      args[n_args++] = level_value;
 
      // // set drop tol for ILU(k) factorization
      // args[n_args++] = "-sparseA";
      // char droptol[20];
      // sprintf(droptol,"%f",drop_tol);
      // args[n_args++] = droptol;
 
      // // set ILUT factorization if required
      // if (use_ilut)
      //  {
      //   args[n_args++] = "-ilut";
      //   args[n_args++] = droptol;
      //  }
 
      // set printing of Euclid data
      if (print_level == 0)
      {
        args[n_args++] = "-eu_stats";
        args[n_args++] = "0";
        args[n_args++] = "-eu_mem";
        args[n_args++] = "0";
      }
      if (print_level == 1)
      {
        args[n_args++] = "-eu_stats";
        args[n_args++] = "1";
        args[n_args++] = "-eu_mem";
        args[n_args++] = "0";
      }
      if (print_level == 2)
      {
        args[n_args++] = "-eu_stats";
        args[n_args++] = "1";
        args[n_args++] = "-eu_mem";
        args[n_args++] = "1";
      }
 
      // set next entry in array to null
      args[n_args] = 0;
 
      // Send the parameters
      HYPRE_EuclidSetParams(euclid_object, n_args, const_cast<char**>(args));
    }

References compute_granudrum_aor::args.

Referenced by oomph::HypreInterface::hypre_solver_setup().

Variable Documentation

AMG_coarsening

unsigned oomph::HypreHelpers::AMG_coarsening = 6

Default AMG coarsening strategy. Coarsening types include: 0 = CLJP (parallel coarsening using independent sets) 1 = classical RS with no boundary treatment (not recommended in parallel) 3 = modified RS with 3rd pass to add C points on the boundaries 6 = Falgout (uses 1 then CLJP using interior coarse points as first independent set) 8 = PMIS (parallel coarsening using independent sets - lower complexities than 0, maybe also slower convergence) 10= HMIS (one pass RS on each processor then PMIS on interior coarse points as first independent set) 11= One pass RS on each processor (not recommended)

Referenced by oomph::HypreInterface::HypreInterface().

AMG_strength

double oomph::HypreHelpers::AMG_strength = 0.25

Default for AMG strength (0.25 recommended for 2D problems; larger (0.5-0.75, say) for 3D

Referenced by oomph::HypreInterface::HypreInterface().

AMG_truncation

double oomph::HypreHelpers::AMG_truncation = 0.0

AMG interpolation truncation factor.

Referenced by oomph::HypreInterface::HypreInterface().