oomph::HSL_MA42 Class Reference

#include <frontal_solver.h>

Inheritance diagram for oomph::HSL_MA42:

Public Member Functions
	HSL_MA42 ()
	~HSL_MA42 ()
	Destructor, clean up the allocated memory. More...
	HSL_MA42 (const HSL_MA42 &)=delete
	Broken copy constructor. More...
void	operator= (const HSL_MA42 &)=delete
	Broken assignment operator. More...
void	clean_up_memory ()
	Clean up memory. More...
void	disable_resolve ()
	Overload disable resolve so that it cleans up memory too. More...
void	solve (Problem *const &problem_pt, DoubleVector &result)
	Wrapper for HSL MA42 frontal solver. More...
void	solve (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
void	solve (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
void	resolve (const DoubleVector &rhs, DoubleVector &result)
	Wrapper for HSL MA42 frontal solver. More...
void	reorder_elements (Problem *const &problem_pt)
	Function to reorder the elements based on Sloan's algorithm. More...
void	enable_doc_stats ()
	Enable documentation of statistics. More...
void	disable_doc_stats ()
	Disable documentation of statistics. More...
void	enable_reordering ()
	Enable reordering using Sloan's algorithm. More...
void	disable_reordering ()
	Disable reordering. More...
void	enable_direct_access_files ()
	Enable use of direct access files. More...
void	disable_direct_access_files ()
	Disable use of direct access files. More...
double &	lenbuf_factor0 ()
double &	lenbuf_factor1 ()
double &	lenbuf_factor2 ()
double &	front_factor ()
double &	lenfle_factor ()
Public Member Functions inherited from oomph::LinearSolver
	LinearSolver ()
	Empty constructor, initialise the member data. More...
	LinearSolver (const LinearSolver &dummy)=delete
	Broken copy constructor. More...
void	operator= (const LinearSolver &)=delete
	Broken assignment operator. More...
virtual	~LinearSolver ()
	Empty virtual destructor. More...
void	enable_doc_time ()
	Enable documentation of solve times. More...
void	disable_doc_time ()
	Disable documentation of solve times. More...
bool	is_doc_time_enabled () const
	Is documentation of solve times enabled? More...
bool	is_resolve_enabled () const
	Boolean flag indicating if resolves are enabled. More...
virtual void	enable_resolve ()
virtual void	solve_transpose (Problem *const &problem_pt, DoubleVector &result)
virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
virtual void	resolve_transpose (const DoubleVector &rhs, DoubleVector &result)
virtual double	jacobian_setup_time () const
virtual double	linear_solver_solution_time () const
virtual void	enable_computation_of_gradient ()
void	disable_computation_of_gradient ()
void	reset_gradient ()
void	get_gradient (DoubleVector &gradient)
	function to access the gradient, provided it has been computed More...
Public Member Functions inherited from oomph::DistributableLinearAlgebraObject
	DistributableLinearAlgebraObject ()
	Default constructor - create a distribution. More...
	DistributableLinearAlgebraObject (const DistributableLinearAlgebraObject &matrix)=delete
	Broken copy constructor. More...
void	operator= (const DistributableLinearAlgebraObject &)=delete
	Broken assignment operator. More...
virtual	~DistributableLinearAlgebraObject ()
	Destructor. More...
LinearAlgebraDistribution *	distribution_pt () const
	access to the LinearAlgebraDistribution More...
unsigned	nrow () const
	access function to the number of global rows. More...
unsigned	nrow_local () const
	access function for the num of local rows on this processor. More...
unsigned	nrow_local (const unsigned &p) const
	access function for the num of local rows on this processor. More...
unsigned	first_row () const
	access function for the first row on this processor More...
unsigned	first_row (const unsigned &p) const
	access function for the first row on this processor More...
bool	distributed () const
	distribution is serial or distributed More...
bool	distribution_built () const
void	build_distribution (const LinearAlgebraDistribution *const dist_pt)
void	build_distribution (const LinearAlgebraDistribution &dist)

Private Member Functions
void	solve_for_one_dof (Problem *const &problem_pt, DoubleVector &result)

Private Attributes
bool	Doc_stats
	Doc the solver stats or stay quiet? More...
bool	Reorder_flag
	Reorder elements with Sloan's algorithm? More...
bool	Use_direct_access_files
	Use direct access files? More...
double	Lenbuf_factor0
double	Lenbuf_factor1
double	Lenbuf_factor2
double	Front_factor
double	Lenfle_factor
int	Icntl [8]
	Control flag for MA42; see MA42 documentation for details. More...
int	Isave [45]
	Control flag for MA42; see MA42 documentation for details. More...
int	Info [23]
	Control flag for MA42; see MA42 documentation for details. More...
double *	W
	Workspace storage for MA42. More...
int	Lw
	Size of the workspace array, W. More...
int *	IW
	Integer workspace storage for MA42. More...
int	Liw
	Size of the integer workspace array. More...
unsigned long	N_dof
	Size of the linear system. More...

Additional Inherited Members
Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void	clear_distribution ()
Protected Attributes inherited from oomph::LinearSolver
bool	Enable_resolve
bool	Doc_time
	Boolean flag that indicates whether the time taken. More...
bool	Compute_gradient
bool	Gradient_has_been_computed
	flag that indicates whether the gradient was computed or not More...
DoubleVector	Gradient_for_glob_conv_newton_solve

Detailed Description

Linear solver class that provides a wrapper to the frontal solver MA42 from the HSL library; see http://www.hsl.rl.ac.uk/.

Constructor & Destructor Documentation

HSL_MA42() [1/2]

oomph::HSL_MA42::HSL_MA42 ( )

inline

Constructor: By default suppress verbose output (stats), don't reorder elements and don't use direct access files

               : W(0), Lw(0), IW(0), Liw(0), N_dof(0)
    {
      Doc_stats = false;
      Reorder_flag = false;
      Use_direct_access_files = false;
 
      // Default values for memory allocation
      Lenbuf_factor0 = 1.2;
      Lenbuf_factor1 = 1.2;
      Lenbuf_factor2 = 1.5;
      Front_factor = 1.1;
      Lenfle_factor = 1.5;
    }

References Doc_stats, Front_factor, Lenbuf_factor0, Lenbuf_factor1, Lenbuf_factor2, Lenfle_factor, Reorder_flag, and Use_direct_access_files.

~HSL_MA42()

oomph::HSL_MA42::~HSL_MA42 ( )

inline

Destructor, clean up the allocated memory.

    {
      clean_up_memory();
    }

References clean_up_memory().

HSL_MA42() [2/2]

oomph::HSL_MA42::HSL_MA42 ( const HSL_MA42 &  )

delete

Broken copy constructor.

Member Function Documentation

clean_up_memory()

void oomph::HSL_MA42::clean_up_memory ( )

inlinevirtual

Clean up memory.

Reimplemented from oomph::LinearSolver.

    {
      if (IW)
      {
        delete[] IW;
        IW = 0;
        Liw = 0;
      }
      if (W)
      {
        delete[] W;
        W = 0;
        Lw = 0;
      }
    }

References IW, Liw, Lw, and W.

Referenced by disable_resolve(), and ~HSL_MA42().

disable_direct_access_files()

void oomph::HSL_MA42::disable_direct_access_files ( )

inline

Disable use of direct access files.

    {
      Use_direct_access_files = false;
    }

References Use_direct_access_files.

disable_doc_stats()

void oomph::HSL_MA42::disable_doc_stats ( )

inline

Disable documentation of statistics.

    {
      Doc_stats = false;
    }

References Doc_stats.

disable_reordering()

void oomph::HSL_MA42::disable_reordering ( )

inline

Disable reordering.

    {
      Reorder_flag = false;
    }

References Reorder_flag.

disable_resolve()

void oomph::HSL_MA42::disable_resolve ( )

inlinevirtual

Overload disable resolve so that it cleans up memory too.

Reimplemented from oomph::LinearSolver.

    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }

References clean_up_memory(), and oomph::LinearSolver::disable_resolve().

enable_direct_access_files()

void oomph::HSL_MA42::enable_direct_access_files ( )

inline

Enable use of direct access files.

    {
      Use_direct_access_files = true;
    }

References Use_direct_access_files.

enable_doc_stats()

void oomph::HSL_MA42::enable_doc_stats ( )

inline

Enable documentation of statistics.

    {
      Doc_stats = true;
    }

References Doc_stats.

enable_reordering()

void oomph::HSL_MA42::enable_reordering ( )

inline

Enable reordering using Sloan's algorithm.

    {
      Reorder_flag = true;
    }

References Reorder_flag.

front_factor()

double& oomph::HSL_MA42::front_factor ( )

inline

Factor to increase storage for front size; see MA42 documentation for details.

    {
      return Front_factor;
    }

References Front_factor.

lenbuf_factor0()

double& oomph::HSL_MA42::lenbuf_factor0 ( )

inline

Factor to increase storage for lenbuf[0]; see MA42 documentation for details.

    {
      return Lenbuf_factor0;
    }

References Lenbuf_factor0.

lenbuf_factor1()

double& oomph::HSL_MA42::lenbuf_factor1 ( )

inline

Factor to increase storage for lenbuf[1]; see MA42 documentation for details.

    {
      return Lenbuf_factor1;
    }

References Lenbuf_factor1.

lenbuf_factor2()

double& oomph::HSL_MA42::lenbuf_factor2 ( )

inline

Factor to increase storage for lenbuf[2]; see MA42 documentation for details.

    {
      return Lenbuf_factor2;
    }

References Lenbuf_factor2.

lenfle_factor()

double& oomph::HSL_MA42::lenfle_factor ( )

inline

Factor to increase the size of the direct access files; see MA42 documentation for details.

    {
      return Lenfle_factor;
    }

References Lenfle_factor.

operator=()

void oomph::HSL_MA42::operator= ( const HSL_MA42 &  )

delete

Broken assignment operator.

reorder_elements()

void oomph::HSL_MA42::reorder_elements

(

Problem *const &

problem_pt

)

Function to reorder the elements based on Sloan's algorithm.

Function to reorder the elements according to Sloan's algorithm.

  {
    // Find the number of elements
    unsigned long n_el = problem_pt->mesh_pt()->nelement();
 
    // Find the number of dofs
    int n_dof = problem_pt->ndof();
 
    Vector<int> order(n_el);
 
    // Control parameters
    int icntl[10], info[15], direct, nsup, vars[n_dof], svar[n_dof];
    int liw, lw, *perm = 0;
    double wt[3], rinfo[6];
 
 
    // Direct ordering: 1
    direct = 1;
 
    // Initial guess for workspaces (deliberately too small so the
    // routine fails and returns the required sizes
    liw = 1;
    lw = 1;
 
    // Initialise things here
    MC63ID(icntl);
 
 
    // Suppress printing of error and warning messages?
    if (!Doc_stats)
    {
      icntl[0] = -1;
      icntl[1] = -1;
    }
 
    // Set up iw and w
    int* iw = new int[liw];
    double* w = new double[lw];
 
    // Set up the vectors that hold the element connectivity information
    // Can initialise the first entry of eltptr to 1
    Vector<int> eltvar, eltptr(1, 1);
 
    // Locally cache pointer to assembly handler
    AssemblyHandler* const assembly_handler_pt =
      problem_pt->assembly_handler_pt();
 
    // Now loop over all the elements and assemble eltvar and eltptr
    for (unsigned long e = 0; e < n_el; e++)
    {
      // Set up the default order
      order[e] = e;
 
      // Get pointer to the element
      GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
 
      // Get the number of variables in the element via the assembly handler
      int nvar = assembly_handler_pt->ndof(elem_pt);
 
      // Multiple equations
      //-------------------
      if (nvar != 1)
      {
        // Copy the equation numbers into the global array
        for (int i = 0; i < nvar; i++)
        {
          // Need to add one to equation numbers
          eltvar.push_back(1 + assembly_handler_pt->eqn_number(elem_pt, i));
        }
        eltptr.push_back(eltptr.back() + nvar);
      }
      // Just one equation
      //------------------
      else
      {
        // Here's the number of the only equation
        int only_eqn = assembly_handler_pt->eqn_number(elem_pt, 0);
 
        // Need to add one to equation numbers
        eltvar.push_back(1 + only_eqn);
 
        // Now add a dummy eqn either before or after the current one
        int dummy_eqn;
        if (only_eqn != (n_dof - 1))
        {
          dummy_eqn = only_eqn + 1;
        }
        else
        {
          dummy_eqn = only_eqn - 1;
        }
 
        eltvar.push_back(1 + dummy_eqn);
 
        eltptr.push_back(eltptr.back() + 2);
      }
 
    } // End of loop over the elements
 
    // Set the value of n_e (number of entries in the element variable list
    int n_e = eltvar.size();
 
    do
    {
      // Call the reordering routine
      MC63AD(direct,
             n_dof,
             n_el,
             n_e,
             &eltvar[0],
             &eltptr[0],
             &order[0],
             perm,
             nsup,
             vars,
             svar,
             wt,
             liw,
             iw,
             lw,
             w,
             icntl,
             info,
             rinfo);
 
      // If not enough space in iw or w, allocate enough and try again
      if (info[0] == -4)
      {
        if (Doc_stats)
          oomph_info << " Reallocating liw to " << info[4] << std::endl;
        delete[] iw;
        liw = info[4];
        iw = new int[liw];
      }
 
      // If not enough space in w, allocate enough and try again
      if (info[0] == -2)
      {
        if (Doc_stats)
          oomph_info << " Reallocating lw to " << info[5] << std::endl;
        delete[] w;
        lw = info[5];
        w = new double[lw];
      }
 
    } while (info[0] < 0);
 
 
    if (Doc_stats)
    {
      oomph_info << "\n Initial wavefront details  :\n      max " << rinfo[0]
                 << "\tmean " << rinfo[1] << "\tprofile " << rinfo[2];
      oomph_info << "\n Reordered wavefront details:\n      max " << rinfo[3]
                 << "\tmean " << rinfo[4] << "\tprofile " << rinfo[5];
      oomph_info << std::endl;
    }
 
    // Free the memory allocated
    delete[] w;
    delete[] iw;
 
 
    // Store re-ordered elements in temporary vector
    Vector<GeneralisedElement*> tmp_element_pt(n_el);
    for (unsigned e = 0; e < n_el; e++)
    {
      tmp_element_pt[e] = problem_pt->mesh_pt()->element_pt(order[e] - 1);
    }
    for (unsigned e = 0; e < n_el; e++)
    {
      problem_pt->mesh_pt()->element_pt(e) = tmp_element_pt[e];
    }
  }

References oomph::Problem::assembly_handler_pt(), Doc_stats, e(), oomph::Mesh::element_pt(), oomph::AssemblyHandler::eqn_number(), i, info, MC63AD, MC63ID, oomph::Problem::mesh_pt(), oomph::Problem::ndof(), oomph::AssemblyHandler::ndof(), oomph::Mesh::nelement(), oomph::oomph_info, and w.

Referenced by solve().

resolve()

void oomph::HSL_MA42::resolve ( const DoubleVector &  rhs, DoubleVector &  result  )

virtual

Wrapper for HSL MA42 frontal solver.

Return the solution to the linear system Ax = result, where A is the most recently factorised jacobian matrix of the problem problem_pt. The solution is returned in the result vector.

Reimplemented from oomph::LinearSolver.

  {
    // If we haven't stored the factors complain
    if (W == 0)
    {
      throw OomphLibError("Resolve called without first solving",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // Find the number of dofs (variables)
    int n_dof = N_dof;
 
    // Check that the number of DOFs is equal to the size of the RHS vector
#ifdef PARANOID
    if (n_dof != static_cast<int>(rhs.nrow()))
    {
      throw OomphLibError(
        "RHS does not have the same dimension as the linear system",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Special case: just one variable! MA42 can't handle it
    // Solve using the stored jacobian (single double)
    if (n_dof == 1)
    {
      result[0] = rhs[0] / W[0];
      return;
    }
 
    // Otherwise actually call the MA42 routine
    // We are solving the original matrix (not the transpose)
    int trans = 0;
    // There is only one rhs
    int nrhs = 1;
    // The size of the vectors is the number of degrees of freedom in the
    // problem
    int lx = n_dof;
 
    // Allocate storage for the RHS
    double** b = new double*;
    *b = new double[n_dof];
    // Load in the RHS vector
    for (int i = 0; i < n_dof; i++)
    {
      b[0][i] = rhs[i];
    }
 
    // Storage for the results
    double** x = new double*;
    *x = new double[n_dof];
 
    // Now call the resolver
    MA42CD(trans, nrhs, lx, b, x, Lw, W, Liw, IW, Icntl, Isave, Info);
 
    // If there has been an error throw it
    if (Info[0] < 0)
    {
      throw NewtonSolverError(true);
    }
 
    result.build(this->distribution_pt());
    // Put the answer into the result array
    for (int i = 0; i < n_dof; ++i)
    {
      result[i] = x[0][i];
    }
 
    // Delete the allocated storage
    delete[] * x;
    delete x;
    delete[] * b;
    delete b;
  }

References b, oomph::DoubleVector::build(), oomph::DistributableLinearAlgebraObject::distribution_pt(), i, Icntl, Info, Isave, IW, Liw, Lw, Mesh_Parameters::lx, MA42CD, N_dof, oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, trans, W, and plotDoE::x.

solve() [1/3]

void oomph::HSL_MA42::solve ( DoubleMatrixBase *const &  matrix_pt, const DoubleVector &  rhs, DoubleVector &  result  )

inlinevirtual

Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the linear system. Call the broken base-class version. If you want this, please implement it

Reimplemented from oomph::LinearSolver.

    {
      LinearSolver::solve(matrix_pt, rhs, result);
    }

References oomph::LinearSolver::solve().

solve() [2/3]

void oomph::HSL_MA42::solve ( DoubleMatrixBase *const &  matrix_pt, const Vector< double > &  rhs, Vector< double > &  result  )

inlinevirtual

Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the linear system Call the broken base-class version. If you want this, please implement it

Reimplemented from oomph::LinearSolver.

    {
      LinearSolver::solve(matrix_pt, rhs, result);
    }

References oomph::LinearSolver::solve().

solve() [3/3]

void oomph::HSL_MA42::solve ( Problem *const &  problem_pt, DoubleVector &  result  )

virtual

Wrapper for HSL MA42 frontal solver.

Solver: Takes pointer to problem and returns the results Vector which contains the solution of the linear system defined by the problem's fully assembled Jacobian and residual Vector.

Implements oomph::LinearSolver.

  {
#ifdef OOMPH_HAS_MPI
    if (problem_pt->communicator_pt()->nproc() > 1)
    {
      std::ostringstream error_message;
      error_message
        << "HSL_MA42 solver cannot be used in parallel.\n"
        << "Please change to another linear solver.\n"
        << "If you want to use a frontal solver then try MumpsSolver\n";
 
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
 
    // Find the number of elements
    unsigned long n_el = problem_pt->mesh_pt()->nelement();
 
    // Find the number of dofs (variables) and store for resolves
    N_dof = problem_pt->ndof();
 
    // Build the distribution .. this is a serial solver
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), N_dof, false);
    this->build_distribution(dist);
 
    // Cast the number of dofs into an integer for the HSL solver
    int n_dof = N_dof;
 
    // Special case: Just one variable! MA42 can't handle that
    if (n_dof == 1)
    {
      solve_for_one_dof(problem_pt, result);
      return;
    }
 
    // Reorder?
    if (this->Reorder_flag)
    {
      reorder_elements(problem_pt);
    }
 
    // Set the control flags
    int ifsize[5];
    double cntl[2], rinfo[2];
 
    // Set up memory for last and for ndf
    int ndf, nmaxe = 2, nrhs = 1, lrhs = 1;
 
    // Memory for last and dx should be allocated dynamically from the heap
    // rather than the stack, otherwise one gets a nasty segmentation fault
    int* last = new int[n_dof];
    // Set up memory for dx
    double** dx = new double*;
    *dx = new double[n_dof];
 
    // Provide storage for exact values required for lenbuf array
    int exact_lenbuf[3] = {0, 0, 0};
    bool exact_lenbuf_available = false;
 
    // Storage for recommendations for lenbuf so we can check how
    // close our factors are
    int lenbuf0_recommendation = 0;
    int lenbuf1_recommendation = 0;
    int lenbuf2_recommendation = 0;
 
    // Provide storage
    int lenbuf[3];
 
 
    // Provide storage for exact values required for lenfle array
    int exact_lenfle[3] = {0, 0, 0};
    bool exact_lenfle_available = false;
 
    // Storage for recommendation for lenfle so we can check how
    // close our factor is
    int lenfle0_recommendation = 0;
    int lenfle1_recommendation = 0;
    int lenfle2_recommendation = 0;
 
    // Provide storage
    int lenfle[3];
 
 
    // Storage for recommendations for front size so we can check how
    // close our factors are
    double front0_recommendation = 0;
    double front1_recommendation = 0;
 
 
    // Flag to indicate if exact, required values for nfront are available
    // from previous unsucessful solve
    bool exact_nfront_available = false;
 
    // Storage for front sizes
    int nfront[2];
 
    // Locally cache pointer to assembly handler
    AssemblyHandler* const assembly_handler_pt =
      problem_pt->assembly_handler_pt();
 
 
    // Loop over size allocation/solver until we've made all the arrays
    //=================================================================
    // big enough...
    //==============
    bool not_done = true;
    while (not_done)
    {
      // Initialise frontal solver
      //=========================
      MA42ID(Icntl, cntl, Isave);
 
 
      // Loop over the elements to setup variables associated with elements
      //==================================================================
      for (unsigned long e = 0; e < n_el; e++)
      {
        // Pointer to the element
        GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
 
        // MH: changed array to allocateable
        int* ivar;
 
        // Get number of variables in element
        int nvar = assembly_handler_pt->ndof(elem_pt);
 
 
        // Multiple equations
        //-------------------
        if (nvar != 1)
        {
          // Copy global equation numbers into local array
          ivar = new int[nvar];
 
          for (int i = 0; i < nvar; i++)
          {
            // Need to add one to equation numbers
            ivar[i] = 1 + assembly_handler_pt->eqn_number(elem_pt, i);
          }
        }
 
        // Just one equation
        //------------------
        else
        {
          // Copy global equation numbers into local array - minimum length: 2
          ivar = new int[2];
 
          // Here's the number of the only equation
          int only_eqn = assembly_handler_pt->eqn_number(elem_pt, 0);
 
          // Need to add one to equation numbers
          ivar[0] = 1 + only_eqn;
 
          // Now add a dummy eqn either before or after the current one
          int dummy_eqn;
          if (only_eqn != (n_dof - 1))
          {
            dummy_eqn = only_eqn + 1;
          }
          else
          {
            dummy_eqn = only_eqn - 1;
          }
          ivar[1] = 1 + dummy_eqn;
 
          nvar = 2;
        }
 
 
        // Now call the frontal routine
        // GB: added check to exclude case with no dofs
        if (nvar)
        {
          MA42AD(nvar, ivar, ndf, last, n_dof, Icntl, Isave, Info);
          // ALH : throw an exception if the frontal_solver fails
          if (Info[0] < 0)
          {
            throw NewtonSolverError(true);
          }
        }
 
        // Cleanup
        delete[] ivar;
        ivar = 0;
      }
 
 
      // Loop over the elements again for symbolic factorisation
      //=======================================================
      for (unsigned long e = 0; e < n_el; e++)
      {
        // Pointer to the element
        GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
 
        // MH: changed array to allocateable
        int* ivar;
 
        // Get number of variables in element
        int nvar = assembly_handler_pt->ndof(elem_pt);
 
        // Multiple equations
        //-------------------
        if (nvar != 1)
        {
          // Copy global equation numbers into local array
          ivar = new int[nvar];
 
          for (int i = 0; i < nvar; i++)
          {
            // Need to add one to equation numbers
            ivar[i] = 1 + assembly_handler_pt->eqn_number(elem_pt, i);
          }
        }
 
        // Just one equation
        //------------------
        else
        {
          // Copy global equation numbers into local array - minimum length: 2
          ivar = new int[2];
 
          // Here's the number of the only equation
          int only_eqn = assembly_handler_pt->eqn_number(elem_pt, 0);
 
          // Need to add one to equation numbers
          ivar[0] = 1 + only_eqn;
 
          // Now add a dummy eqn either before or after the current one
          int dummy_eqn;
          if (only_eqn != (n_dof - 1))
          {
            dummy_eqn = only_eqn + 1;
          }
          else
          {
            dummy_eqn = only_eqn - 1;
          }
          ivar[1] = 1 + dummy_eqn;
 
          nvar = 2;
        }
 
        // GB: added check to exclude case with no dofs
        if (nvar)
        {
          MA42JD(nvar, ivar, ndf, last, nmaxe, ifsize, Icntl, Isave, Info);
 
          // ALH : throw an exception if the frontal_solver fails
          if (Info[0] < 0)
          {
            throw NewtonSolverError(true);
          }
        }
 
        // Cleanup
        delete[] ivar;
        ivar = 0;
 
      } // end of symbolic factorisation
 
 
      // Front sizes: "Somewhat larger than..."
      //---------------------------------------
      front0_recommendation = ifsize[0];
      front1_recommendation = ifsize[1];
      if (!exact_nfront_available)
      {
        nfront[0] = int(ceil(Front_factor * double(front0_recommendation)));
        nfront[1] = int(ceil(Front_factor * double(front1_recommendation)));
 
        if (n_dof < nfront[0]) nfront[0] = n_dof;
        if (n_dof < nfront[1]) nfront[1] = n_dof;
      }
      // We have a problem if the front size is cocked
 
 
      // Sizes for lenbuf[]: "Somewhat larger than..."
      //----------------------------------------------
      lenbuf0_recommendation = ifsize[2] + ndf;
      // If we are storing the matrix,
      // we need to allocate storage for lenbuf_1
      if (Enable_resolve)
      {
        lenbuf1_recommendation = ifsize[3];
      }
      // Otherwise set it to zero
      else
      {
        lenbuf1_recommendation = 0;
      }
      lenbuf2_recommendation = ifsize[4];
 
 
      // Enable use of direct access files?
      //----------------------------------
      if (Use_direct_access_files)
      {
        if (Doc_stats)
        {
          oomph_info << "Using direct access files" << std::endl;
        }
 
        // Unit numbers for direct access files (middle one only used for
        // re-solve; set to zero as this is not used here).
        int istrm[3];
        istrm[0] = 17;
        // If we are storing the matrix, set the stream
        if (Enable_resolve)
        {
          istrm[1] = 19;
        }
        // otherwise, set it to zero
        else
        {
          istrm[1] = 0;
        }
        istrm[2] = 18;
 
        // Set up the memory: Need memory "somewhat larger" than ...
        double factor = 1.1;
        lenbuf[0] = int(ceil(factor * double(10 * (nfront[0] + 1))));
        // If we are storing the matrix, set the value
        if (Enable_resolve)
        {
          lenbuf[1] = int(ceil(factor * double(10 * (nfront[0]))));
        }
        // Otherwise, set to zero
        else
        {
          lenbuf[1] = 0;
        }
        lenbuf[2] = int(ceil(factor * double(10 * (2 * nfront[0] + 5))));
 
 
        // Initial size allocation: Need memory "somewhat larger" than ...
        if (exact_lenfle_available)
        {
          lenfle[0] = exact_lenfle[0];
          lenfle[1] = exact_lenfle[1];
          lenfle[2] = exact_lenfle[2];
        }
        else
        {
          lenfle0_recommendation = ifsize[2] + ndf;
          lenfle1_recommendation = ifsize[3];
          lenfle2_recommendation = ifsize[4];
          lenfle[0] = int(ceil(Lenfle_factor * double(lenfle0_recommendation)));
          lenfle[1] = int(ceil(Lenfle_factor * double(lenfle1_recommendation)));
          lenfle[2] = int(ceil(Lenfle_factor * double(lenfle2_recommendation)));
        }
 
        // Setup direct access files
        MA42PD(istrm, lenbuf, lenfle, Icntl, Isave, Info);
      }
      else
      {
        if (Doc_stats)
        {
          oomph_info << "Not using direct access files" << std::endl;
        }
 
        // Initial size allocation: Need memory "somewhat larger" than ...
        if (exact_lenbuf_available)
        {
          lenbuf[0] = exact_lenbuf[0];
          lenbuf[1] = exact_lenbuf[1];
          lenbuf[2] = exact_lenbuf[2];
        }
        else
        {
          lenbuf[0] =
            int(ceil(Lenbuf_factor0 * double(lenbuf0_recommendation)));
          // If we are storing the matrix, set the buffer size
          if (Enable_resolve)
          {
            lenbuf[1] =
              int(ceil(Lenbuf_factor1 * double(lenbuf1_recommendation)));
          }
          // Otherwise its zero
          else
          {
            lenbuf[1] = 0;
          }
          lenbuf[2] =
            int(ceil(Lenbuf_factor2 * double(lenbuf2_recommendation)));
        }
      }
 
 
      if (Doc_stats)
      {
        oomph_info << "\n FRONT SIZE (min and actual): " << ifsize[0] << " "
                   << nfront[0] << std::endl;
      }
 
 
      // Initialise success
      bool success = true;
 
      // Workspace arrays: calculate amount in local variables initiall
      int lw = 1 + lenbuf[0] + lenbuf[1] + nfront[0] * nfront[1];
      if (lrhs * nfront[0] > nrhs * nfront[1])
      {
        lw += lrhs * nfront[0];
      }
      else
      {
        lw += nrhs * nfront[1];
      }
 
      int liw = lenbuf[2] + 2 * nfront[0] + 4 * nfront[1];
 
      // Set up memory for w and iw
      // Again allocate dynamically from the heap, rather than the stack
      // If the required amount of storage has changed, reallocate
      // and store the values in the object member data
      if (lw != Lw)
      {
        // Set the new value of Lw (member data)
        Lw = lw;
        // Delete the previously allocated storage
        if (W) delete[] W;
        // Now allocate new storage
        W = new (std::nothrow) double[Lw];
        if (!W)
        {
          throw OomphLibError("Out of memory in allocation of w\n",
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // If the required amount of storage has changed, reallocate
      if (liw != Liw)
      {
        // Set the new value of Liw (member data)
        Liw = liw;
        // Delete the previously allocated storgae
        if (IW != 0)
        {
          delete[] IW;
        }
        // Now allocate new storage
        IW = new (std::nothrow) int[Liw];
        if (!IW)
        {
          throw OomphLibError("Out of memory in allocation of iw\n",
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Doc
      if (Doc_stats)
      {
        double temp = (lenbuf[0] * sizeof(double) + lenbuf[2] * sizeof(int)) /
                      (1024.0 * 1024.0);
        oomph_info << "\n ESTIMATED MEMORY USAGE " << temp << "Mb" << std::endl;
      }
 
 
      // Now do the actual frontal assembly and solve loop
      //=================================================
      for (unsigned long e = 0; e < n_el; e++)
      {
        // Get pointer to the element
        GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
 
        // Get number of variables in element
        int nvar = assembly_handler_pt->ndof(elem_pt);
 
        // MH: changed array to allocatable
        int* ivar;
 
        // Multiple equations
        //-------------------
        if (nvar != 1)
        {
          // Copy global equation numbers into local array
          ivar = new int[nvar];
 
          for (int i = 0; i < nvar; i++)
          {
            // Need to add one to equation numbers
            ivar[i] = 1 + assembly_handler_pt->eqn_number(elem_pt, i);
          }
          nmaxe = nvar;
        }
 
        // Just one equation
        //------------------
        else
        {
          // Copy global equation numbers into local array - minimum length: 2
          ivar = new int[2];
 
          // Here's the number of the only equation
          int only_eqn = assembly_handler_pt->eqn_number(elem_pt, 0);
 
          // Need to add one to equation numbers
          ivar[0] = 1 + only_eqn;
 
          // Now add a dummy eqn either before or after the current one
          int dummy_eqn;
          if (only_eqn != (n_dof - 1))
          {
            dummy_eqn = only_eqn + 1;
          }
          else
          {
            dummy_eqn = only_eqn - 1;
          }
          ivar[1] = 1 + dummy_eqn;
 
          nmaxe = 2;
        }
 
        // Set up matrices and Vector for call to get_residuals
        Vector<double> residuals(nvar);
        DenseMatrix<double> jacobian(nvar);
 
        // Get the residuals and jacobian
        assembly_handler_pt->get_jacobian(elem_pt, residuals, jacobian);
 
        // Add dummy rows and columns if required
        if (nvar == 1)
        {
          double onlyjac = jacobian(0, 0);
          jacobian.resize(2, 2);
          jacobian(0, 0) = onlyjac;
          jacobian(1, 0) = 0.0;
          jacobian(0, 1) = 0.0;
          jacobian(1, 1) = 0.0;
 
          double onlyres = residuals[0];
          residuals.resize(2);
          residuals[0] = onlyres;
          residuals[1] = 0.0;
 
          nvar = 2;
        }
 
 
        // GB: added check to exclude case with no dofs
        if (nvar)
        {
          // Now translate the residuals and jacobian into form to be
          // taken by stupid FORTRAN frontal solver
          // double avar[nvar][nmaxe], rhs[1][nmaxe];
          // Assign memory on the heap rather than the stack because the ndofs
          // could get too large
          double** avar = new double*[nvar];
          double* tmp = new double[nvar * nmaxe];
          for (int i = 0; i < nvar; i++)
          {
            avar[i] = &tmp[i * nmaxe];
          }
          double** rhs = new double*[1];
          rhs[0] = new double[nmaxe];
 
 
          for (int i = 0; i < nmaxe; i++)
          {
            rhs[0][i] = residuals[i];
            for (int j = 0; j < nvar; j++)
            {
              // Note need to transpose here
              avar[j][i] = jacobian(i, j);
            }
          }
 
          // Call the frontal solver
          MA42BD(nvar,
                 ivar,
                 ndf,
                 last,
                 nmaxe,
                 avar,
                 nrhs,
                 rhs,
                 lrhs,
                 n_dof,
                 dx,
                 nfront,
                 lenbuf,
                 Lw,
                 W,
                 Liw,
                 IW,
                 Icntl,
                 cntl,
                 Isave,
                 Info,
                 rinfo);
 
 
          // Error check:
          if (Info[0] < 0)
          {
            if (Doc_stats)
              oomph_info << "Error: Info[0]=" << Info[0] << std::endl;
 
            // Solve isn't going to be successful
            success = false;
 
            // Error: Entries in lenbuf too small -- this is covered
            if (Info[0] == -16)
            {
              if (Doc_stats)
                oomph_info << "lenbuf[] too small -- can recover..."
                           << std::endl;
            }
            else if (Info[0] == -12)
            {
              if (Doc_stats)
                oomph_info << "nfront[] too small -- can recover..."
                           << std::endl;
            }
            else if (Info[0] == -17)
            {
              if (Doc_stats)
                oomph_info << "lenfle[] too small -- can recover..."
                           << std::endl;
            }
            else
            {
              if (Doc_stats)
              {
                oomph_info << "Can't recover from this error" << std::endl;
              }
              throw NewtonSolverError(true);
            }
          }
 
          // Clean up the memory
          delete[] rhs[0];
          rhs[0] = 0;
          delete[] rhs;
          rhs = 0;
          delete[] avar[0];
          avar[0] = 0;
          delete[] avar;
        }
 
        // Cleanup
        delete[] ivar;
        ivar = 0;
 
      } // end of loop over elements for assemble/solve
 
      // Set the sign of the jacobian matrix
      problem_pt->sign_of_jacobian() = Info[1];
 
      // If we are not storing the matrix, then delete the assigned workspace
      if (!Enable_resolve)
      {
        // Free the workspace memory assigned and set stored values to zero
        delete[] IW;
        IW = 0;
        Liw = 0;
        delete[] W;
        W = 0;
        Lw = 0;
        // Reset the number of dofs
        N_dof = 0;
      }
 
      // We've done the elements -- did we encounter an error
      // that forces us to re-allocate workspace?
      if (success)
      {
        // We're done!
        not_done = false;
      }
      else
      {
        // Reset sizes for lenbuf
        if (!Use_direct_access_files)
        {
          exact_lenbuf[0] = Info[4];
          exact_lenbuf[1] = Info[5];
          exact_lenbuf[2] = Info[6];
          exact_lenbuf_available = true;
        }
 
        // nfront[] has already been overwritten with required values -- don't
        // need to update anything. Just indicate that new values are
        // available so they don't have to be re-assigned.
        exact_nfront_available = true;
 
        // Reset sizes for lenfle
        exact_lenfle[0] = lenbuf[0] * Info[9];
        exact_lenfle[1] = lenbuf[1] * Info[10];
        exact_lenfle[2] = lenbuf[2] * Info[11];
        exact_lenfle_available = true;
      }
 
    } // end of loop over shuffling of workspace array sizes until it all
      // fits...
 
    result.build(&dist);
    // Now load the values back into result
    for (int i = 0; i < n_dof; i++) result[i] = dx[0][i];
 
    // Print the actual memory used
    if (Doc_stats)
    {
      if (!Use_direct_access_files)
      {
        double temp = (Info[4] * sizeof(double) + Info[6] * sizeof(int)) /
                      (1024.0 * 1024.0);
        oomph_info << " TOTAL MEMORY USED " << temp << "Mb" << std::endl;
      }
      oomph_info << std::endl;
      oomph_info << "lenbuf[]: " << lenbuf[0] << " " << lenbuf[1] << " "
                 << lenbuf[2] << " " << std::endl;
      oomph_info << "lenbuf[] factors required and initially specified:"
                 << std::endl;
      oomph_info << "lenbuf[0]: " << Info[4] / double(lenbuf0_recommendation)
                 << " " << Lenbuf_factor0 << std::endl;
      oomph_info << "lenbuf[1]: " << Info[5] / double(lenbuf1_recommendation)
                 << " " << Lenbuf_factor1 << std::endl;
 
      oomph_info << "lenbuf[2]: " << Info[6] / double(lenbuf2_recommendation)
                 << " " << Lenbuf_factor2 << std::endl;
      oomph_info << "nfront[] factors required and initially specified:"
                 << std::endl;
      oomph_info << "nfront[0]: " << nfront[0] / double(front0_recommendation)
                 << " " << Front_factor << std::endl;
      oomph_info << "nfront[1]: " << nfront[1] / double(front1_recommendation)
                 << " " << Front_factor << std::endl;
      if (Use_direct_access_files)
      {
        oomph_info << "lenfle[]: " << lenfle[0] << " " << lenfle[1] << " "
                   << lenfle[2] << " " << std::endl;
        oomph_info << "lenfle[] factors required and initially specified:"
                   << std::endl;
        oomph_info << "lenfle[0]: "
                   << Info[9] * lenbuf[0] / double(lenfle0_recommendation)
                   << " " << Lenfle_factor << std::endl;
        oomph_info << "lenfle[1]: "
                   << Info[10] * lenbuf[1] / double(lenfle1_recommendation)
                   << " " << Lenfle_factor << std::endl;
        oomph_info << "lenfle[2]: "
                   << Info[11] * lenbuf[2] / double(lenfle2_recommendation)
                   << " " << Lenfle_factor << std::endl;
      }
    }
 
    // Free the memory assigned
    delete[] * dx;
    delete dx;
    delete[] last;
  }

References oomph::Problem::assembly_handler_pt(), oomph::DoubleVector::build(), oomph::DistributableLinearAlgebraObject::build_distribution(), Eigen::bfloat16_impl::ceil(), oomph::Problem::communicator_pt(), Doc_stats, e(), oomph::Mesh::element_pt(), oomph::LinearSolver::Enable_resolve, oomph::AssemblyHandler::eqn_number(), Front_factor, oomph::AssemblyHandler::get_jacobian(), i, Icntl, Info, int(), Isave, IW, j, Eigen::placeholders::last, Lenbuf_factor0, Lenbuf_factor1, Lenbuf_factor2, Lenfle_factor, Liw, Lw, MA42AD, MA42BD, MA42ID, MA42JD, MA42PD, oomph::Problem::mesh_pt(), N_dof, oomph::Problem::ndof(), oomph::AssemblyHandler::ndof(), oomph::Mesh::nelement(), oomph::OomphCommunicator::nproc(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, reorder_elements(), Reorder_flag, oomph::DenseMatrix< T >::resize(), oomph::Problem::sign_of_jacobian(), solve_for_one_dof(), tmp, Use_direct_access_files, and W.

solve_for_one_dof()

void oomph::HSL_MA42::solve_for_one_dof ( Problem *const &  problem_pt, DoubleVector &  result  )

private

Special solver for problems with 1 dof (MA42 can't handle this case so solve() forwards the "solve" to this function.

Special solver for problems with one DOF – HSL_MA42 can't handle that!

  {
    // Find the number of elements
    unsigned long n_el = problem_pt->mesh_pt()->nelement();
 
#ifdef PARANOID
    // Find the number of dofs (variables)
    int n_dof = problem_pt->ndof();
    // Check
    if (n_dof != 1)
    {
      throw OomphLibError(
        "You can only call this if the problem has just one dof!",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Global jac and residuals are scalars!
    double global_jac = 0.0;
    double global_res = 0.0;
 
    // Locally cache pointer to assembly handler
    AssemblyHandler* const assembly_handler_pt =
      problem_pt->assembly_handler_pt();
 
 
    // Do the assembly
    for (unsigned long e = 0; e < n_el; e++)
    {
      // Get pointer to the element
      GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
 
      // Get number of variables in element
      int nvar = assembly_handler_pt->ndof(elem_pt);
 
      // Only assemble if there's something to be assembled
      if (nvar > 0)
      {
        // Set up matrices and Vector for call to get_residuals
        Vector<double> residuals(nvar);
        DenseMatrix<double> jacobian(nvar);
 
        // Get the residuals and jacobian
        assembly_handler_pt->get_jacobian(elem_pt, residuals, jacobian);
 
        // Add contribution
        global_jac += jacobian(0, 0);
        global_res += residuals[0];
      }
    }
 
    // "Solve"
    result[0] = global_res / global_jac;
 
    // If we are storing the matrix, allocate the storage
    if (Enable_resolve)
    {
      // If storage has been allocated delete it
      if (W != 0)
      {
        delete[] W;
      }
      // Now allocate new storage
      W = new double[1];
      // Store the global jacobian
      W[0] = global_jac;
      // Set the number of degrees of freedom
      N_dof = 1;
    }
    // Set the sign of the jacobian
    problem_pt->sign_of_jacobian() =
      static_cast<int>(std::fabs(global_jac) / global_jac);
  }

References oomph::Problem::assembly_handler_pt(), e(), oomph::Mesh::element_pt(), oomph::LinearSolver::Enable_resolve, boost::multiprecision::fabs(), oomph::AssemblyHandler::get_jacobian(), oomph::Problem::mesh_pt(), N_dof, oomph::Problem::ndof(), oomph::AssemblyHandler::ndof(), oomph::Mesh::nelement(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::Problem::sign_of_jacobian(), and W.

Referenced by solve().

Member Data Documentation

Doc_stats

bool oomph::HSL_MA42::Doc_stats

private

Doc the solver stats or stay quiet?

Referenced by disable_doc_stats(), enable_doc_stats(), HSL_MA42(), reorder_elements(), and solve().

Front_factor

double oomph::HSL_MA42::Front_factor

private

Factor to increase storage for front size; see MA42 documentation for details.

Referenced by front_factor(), HSL_MA42(), and solve().

Icntl

int oomph::HSL_MA42::Icntl[8]

private

Control flag for MA42; see MA42 documentation for details.

Referenced by resolve(), and solve().

Info

int oomph::HSL_MA42::Info[23]

private

Control flag for MA42; see MA42 documentation for details.

Referenced by resolve(), and solve().

Isave

int oomph::HSL_MA42::Isave[45]

private

Control flag for MA42; see MA42 documentation for details.

Referenced by resolve(), and solve().

IW

int* oomph::HSL_MA42::IW

private

Integer workspace storage for MA42.

Referenced by clean_up_memory(), resolve(), and solve().

Lenbuf_factor0

double oomph::HSL_MA42::Lenbuf_factor0

private

Factor to increase storage for lenbuf[0]; see MA42 documentation for details.

Referenced by HSL_MA42(), lenbuf_factor0(), and solve().

Lenbuf_factor1

double oomph::HSL_MA42::Lenbuf_factor1

private

Factor to increase storage for lenbuf[1]; see MA42 documentation for details

Referenced by HSL_MA42(), lenbuf_factor1(), and solve().

Lenbuf_factor2

double oomph::HSL_MA42::Lenbuf_factor2

private

Factor to increase storage for lenbuf[2]; see MA42 documentation for details.

Referenced by HSL_MA42(), lenbuf_factor2(), and solve().

Lenfle_factor

double oomph::HSL_MA42::Lenfle_factor

private

Factor to increase size of direct access files; see MA42 documentation for details.

Referenced by HSL_MA42(), lenfle_factor(), and solve().

Liw

int oomph::HSL_MA42::Liw

private

Size of the integer workspace array.

Referenced by clean_up_memory(), resolve(), and solve().

Lw

int oomph::HSL_MA42::Lw

private

Size of the workspace array, W.

Referenced by clean_up_memory(), resolve(), and solve().

N_dof

unsigned long oomph::HSL_MA42::N_dof

private

Size of the linear system.

Referenced by resolve(), solve(), and solve_for_one_dof().

Reorder_flag

bool oomph::HSL_MA42::Reorder_flag

private

Reorder elements with Sloan's algorithm?

Referenced by disable_reordering(), enable_reordering(), HSL_MA42(), and solve().

Use_direct_access_files

bool oomph::HSL_MA42::Use_direct_access_files

private

Use direct access files?

Referenced by disable_direct_access_files(), enable_direct_access_files(), HSL_MA42(), and solve().

W

double* oomph::HSL_MA42::W

private

Workspace storage for MA42.

Referenced by clean_up_memory(), resolve(), solve(), and solve_for_one_dof().

---

The documentation for this class was generated from the following files:

• frontal_solver.h
• frontal_solver.cc