assembly_handler.cc

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
//OOOMPH-LIB includes
#include "assembly_handler.h"
#include "elements.h"
#include "problem.h"
#include "mesh.h"

namespace oomph
{
 ///////////////////////////////////////////////////////////////////////
 // Non-inline functions for the AssemblyHandler class
 //////////////////////////////////////////////////////////////////////

 //===================================================================
 ///Get the number of elemental degrees of freedom. Direct call
 ///to the function in the element.
 //===================================================================
 unsigned AssemblyHandler::ndof(GeneralisedElement* const &elem_pt)
 {return elem_pt->ndof();}

 //===================================================================
 ///Get the vector of dofs in the element elem_pt at time level t
 ///Direct call to the function in the element.
 //===================================================================
 void AssemblyHandler::dof_vector(GeneralisedElement* const &elem_pt,
                                 const unsigned &t, Vector<double> &dof)
 {return elem_pt->dof_vector(t,dof);}

 //===================================================================
 ///Get the vector of pointers to dofs in the element elem_pt
 ///Direct call to the function in the element.
 //===================================================================
 void AssemblyHandler::dof_pt_vector(GeneralisedElement* const &elem_pt,
                                     Vector<double*> &dof_pt)
 {return elem_pt->dof_pt_vector(dof_pt);}

 /// \short Return the t-th level of storage associated with the i-th
 /// (local) dof stored in the problem
 double &AssemblyHandler::local_problem_dof(Problem* const &problem_pt, 
                                            const unsigned &t, 
                                            const unsigned &i)
 {return *(problem_pt->dof_pt(i)+t);}
 


 //==================================================================
 ///Get the global equation number of the local unknown. Direct call
 ///to the function in the element.
 //==================================================================
 unsigned long AssemblyHandler::eqn_number(GeneralisedElement* const &elem_pt,
                                           const unsigned &ieqn_local)
 {return elem_pt->eqn_number(ieqn_local);}

 //==================================================================
 ///Get the residuals by calling the underlying element's residuals
 ///directly.
 //=================================================================
 void AssemblyHandler::get_residuals(GeneralisedElement* const &elem_pt,
                                     Vector<double> &residuals)
 {elem_pt->get_residuals(residuals);}

 //=======================================================================
 /// Calculate the elemental Jacobian matrix "d equation / d variable" by
 /// calling the element's get_jacobian function.
 //======================================================================
 void AssemblyHandler::get_jacobian(GeneralisedElement* const &elem_pt,
                                    Vector<double> &residuals,
                                    DenseMatrix<double> &jacobian)
 {elem_pt->get_jacobian(residuals,jacobian);}

 //=======================================================================
 /// Calculate all desired vectors and matrices that are required by
 /// the element by calling the get_jacobian function
 //=======================================================================
 void AssemblyHandler::get_all_vectors_and_matrices(
  GeneralisedElement* const &elem_pt,
  Vector<Vector<double> >&vec, Vector<DenseMatrix<double> > &matrix)
 {
  get_jacobian(elem_pt,vec[0],matrix[0]);
 }

 //=======================================================================
 /// \short Calculate the derivative of the residuals with respect to
 /// a parameter, by calling the elemental function
 //======================================================================
 void AssemblyHandler::get_dresiduals_dparameter(
  GeneralisedElement* const &elem_pt,
  double* const &parameter_pt,
  Vector<double> &dres_dparam)
 {
  elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);
 }


 //=====================================================================
 /// \short Calculate the derivative of the residuals and jacobian
 /// with respect to a parameter by calling the elemental function
 //========================================================================
 void AssemblyHandler::get_djacobian_dparameter(
  GeneralisedElement* const &elem_pt,
  double* const &parameter_pt,
  Vector<double> &dres_dparam,
  DenseMatrix<double> &djac_dparam)
 {
  elem_pt->get_djacobian_dparameter(parameter_pt,dres_dparam,djac_dparam);
 }

 /// \short Calculate the product of the Hessian (derivative of Jacobian with
 /// respect to all variables) an eigenvector, Y, and
 /// other specified vectors, C
 /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
 /// At the moment the dof pointer is passed in to enable
 /// easy calculation of finite difference default
 void AssemblyHandler::get_hessian_vector_products(
  GeneralisedElement* const &elem_pt,
  Vector<double> const &Y,
  DenseMatrix<double> const &C,
  DenseMatrix<double> &product)
 {
  elem_pt->get_hessian_vector_products(Y,C,product);
 }


 //=======================================================================
 /// Return the eigenfunction(s) associated with the bifurcation that
 /// has been detected in bifurcation tracking problems. Default
 /// Broken implementation
 //=======================================================================
 double* AssemblyHandler::bifurcation_parameter_pt() const
 {
  std::ostringstream error_stream;
  error_stream
   <<
   "There is no bifurcation parameter associated with the current assembly handler.\n"
   <<
   "Eigenfunction are only calculated by the Fold, PitchFork and Hopf Handlers"
   << "\n";

  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
  return 0;
 }


 //=======================================================================
 /// Return the eigenfunction(s) associated with the bifurcation that
 /// has been detected in bifurcation tracking problems. Default
 /// Broken implementation
 //=======================================================================
 void AssemblyHandler::get_eigenfunction(
  Vector<DoubleVector> &eigenfunction)
 {
  std::ostringstream error_stream;
  error_stream
   <<
   "There is no eigenfunction associated with the current assembly handler.\n"
   <<
   "Eigenfunction are only calculated by the Fold, PitchFork and Hopf Handlers"
   << "\n";

  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
 }

 ///==========================================================================
 /// Compute the inner products of the given vector of pairs of 
 /// history values over the element. The values of the index in the pair
 /// may be the same.
 //===========================================================================
 void AssemblyHandler::get_inner_products(GeneralisedElement* const &elem_pt,
                                          Vector<std::pair<unsigned,unsigned> >
                                          const &history_index,
                                          Vector<double> &inner_product)
 {
  elem_pt->get_inner_products(history_index,inner_product);
 }

 //==========================================================================
 /// Compute the vectors that when taken as a dot product with
 /// other history values give the inner product over the element.
 /// In other words if we call get_inner_product_vectors({0,1},output)
 /// output[0] will be a vector such that dofs.output[0] gives the inner
 /// product of current dofs with themselves.
 //==========================================================================
 void AssemblyHandler::get_inner_product_vectors(
  GeneralisedElement* const &elem_pt,  Vector<unsigned> const &history_index,
  Vector<Vector<double> >  &inner_product_vector)
 {
  elem_pt->get_inner_product_vectors(history_index,inner_product_vector);
 }



 ///////////////////////////////////////////////////////////////////////
 // Non-inline functions for the ExplicitTimeStepHandler class
 //////////////////////////////////////////////////////////////////////


 //===================================================================
 ///Get the number of elemental degrees of freedom. Direct call
 ///to the function in the element.
 //===================================================================
 unsigned ExplicitTimeStepHandler::ndof(GeneralisedElement* const &elem_pt)
 {return elem_pt->ndof();}

 //==================================================================
 ///Get the global equation number of the local unknown. Direct call
 ///to the function in the element.
 //==================================================================
 unsigned long ExplicitTimeStepHandler::
 eqn_number(GeneralisedElement* const &elem_pt,const unsigned &ieqn_local)
 {return elem_pt->eqn_number(ieqn_local);}

 //==================================================================
 ///Call the element's residuals
 //=================================================================
 void ExplicitTimeStepHandler::get_residuals(
  GeneralisedElement* const &elem_pt,
  Vector<double> &residuals)
 {
  elem_pt->get_residuals(residuals);
 }

 //=======================================================================
 ///Replace get jacobian with the call to get the mass matrix
 //======================================================================
 void ExplicitTimeStepHandler::get_jacobian(GeneralisedElement* const &elem_pt,
                                        Vector<double> &residuals,
                                        DenseMatrix<double> &jacobian)
 {
  elem_pt->get_mass_matrix(residuals,jacobian);
 }


 //=======================================================================
 /// Calculate all desired vectors and matrices that are required by
 /// the problem  by calling those of the underlying element.
 //=======================================================================
 void ExplicitTimeStepHandler::get_all_vectors_and_matrices(
  GeneralisedElement* const &elem_pt,
  Vector<Vector<double> >&vec, Vector<DenseMatrix<double> > &matrix)
 {
#ifdef PARANOID
  //Check dimension
  if(matrix.size() != 1)
   {
    throw OomphLibError(
     "ExplicitTimeSteps should return one matrix",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
#endif
  //Get the residuals and <mass matrix
  elem_pt->get_mass_matrix(vec[0],matrix[0]);
 }



 ///////////////////////////////////////////////////////////////////////
 // Non-inline functions for the EigenProblemHandler class
 //////////////////////////////////////////////////////////////////////


 //===================================================================
 ///Get the number of elemental degrees of freedom. Direct call
 ///to the function in the element.
 //===================================================================
 unsigned EigenProblemHandler::ndof(GeneralisedElement* const &elem_pt)
 {return elem_pt->ndof();}

 //==================================================================
 ///Get the global equation number of the local unknown. Direct call
 ///to the function in the element.
 //==================================================================
 unsigned long EigenProblemHandler::
 eqn_number(GeneralisedElement* const &elem_pt,const unsigned &ieqn_local)
 {return elem_pt->eqn_number(ieqn_local);}

 //==================================================================
 ///Cannot call get_residuals for an eigenproblem, so throw an error
 //=================================================================
 void EigenProblemHandler::get_residuals(GeneralisedElement* const &elem_pt,
                                     Vector<double> &residuals)
 {
  throw OomphLibError(
   "An eigenproblem does not have a get_residuals function",
   OOMPH_CURRENT_FUNCTION,
   OOMPH_EXCEPTION_LOCATION);
 }

 //=======================================================================
 ///Cannot call get_jacobian for an eigenproblem, so throw an error
 //======================================================================
 void EigenProblemHandler::get_jacobian(GeneralisedElement* const &elem_pt,
                                        Vector<double> &residuals,
                                        DenseMatrix<double> &jacobian)
 {
  throw OomphLibError(
   "An eigenproblem does not have a get_jacobian function",
   OOMPH_CURRENT_FUNCTION,
   OOMPH_EXCEPTION_LOCATION);
 }


 //=======================================================================
 /// Calculate all desired vectors and matrices that are required by
 /// the problem  by calling those of the underlying element.
 //=======================================================================
 void EigenProblemHandler::get_all_vectors_and_matrices(
  GeneralisedElement* const &elem_pt,
  Vector<Vector<double> >&vec, Vector<DenseMatrix<double> > &matrix)
 {
#ifdef PARANOID
  //Check dimension
  if(matrix.size() != 2)
   {
    throw OomphLibError("EigenProblems should return two matrices",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  //Find the number of variables
  unsigned n_var = elem_pt->ndof();
  //Assign a dummy residuals vector
  Vector<double> dummy(n_var);
  //Get the jacobian and mass matrices
  elem_pt->get_jacobian_and_mass_matrix(dummy,matrix[0],matrix[1]);

  //If we have a non-zero shift, then shift the A matrix
  if(Sigma_real != 0.0)
   {
    //Set the shifted matrix
    for(unsigned i=0;i<n_var;i++)
     {
      for(unsigned j=0;j<n_var;j++)
       {
        matrix[0](i,j) -= Sigma_real*matrix[1](i,j);
       }
     }
   }
 }


 //======================================================================
 /// Clean up the memory that may have been allocated by the solver
 //=====================================================================
 AugmentedBlockFoldLinearSolver::~AugmentedBlockFoldLinearSolver()
 {
  if(Alpha_pt!=0) {delete Alpha_pt;}
  if(E_pt!=0) {delete E_pt;}
 }

 //===================================================================
 /// Use a block factorisation to solve the augmented system
 /// associated with a PitchFork bifurcation.
 //===================================================================
 void AugmentedBlockFoldLinearSolver::solve(Problem* const &problem_pt,
                                            DoubleVector &result)
 {
  // if the result is setup then it should not be distributed
#ifdef PARANOID
  if (result.built())
   {
    if (result.distributed())
     {
      throw OomphLibError("The result vector must not be distributed",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
#endif

  //Locally cache the pointer to the handler.
  FoldHandler* handler_pt =
   static_cast<FoldHandler*>(problem_pt->assembly_handler_pt());

  //Switch the handler to "block solver" mode
  handler_pt->solve_augmented_block_system();

  //We need to find out the number of dofs in the problem
  unsigned n_dof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
  this->build_distribution(dist);

  // if the result vector is not setup then rebuild with distribution = global
  if (!result.built())
    {
     result.build(this->distribution_pt(),0.0);
    }

  //Setup storage for temporary vectors
  DoubleVector a(this->distribution_pt(),0.0), b(this->distribution_pt(),0.0);

  //Allocate storage for Alpha which can be used in the resolve
  if(Alpha_pt!=0) {delete Alpha_pt;}
  Alpha_pt = new DoubleVector(this->distribution_pt(),0.0);

  //We are going to do resolves using the underlying linear solver
  Linear_solver_pt->enable_resolve();

  //Solve the first system Aa = R
  Linear_solver_pt->solve(problem_pt,a);

  //The vector in the top right-hand block is the jacobian multiplied
  //by the null vector

  //Get the present null vector from the handler
  DoubleVector y(this->distribution_pt(),0.0);
  for(unsigned n=0;n<(n_dof-1);++n)  {y[n] = handler_pt->Y[n];}
  //For simplicity later, add a zero at the end
  y[n_dof-1] = 0.0;

  //Loop over the elements and assemble the product
  //Local storage for the product terms
  DoubleVector Jy(this->distribution_pt(),0.0);

  //Calculate the product of the jacobian matrices, etc
  unsigned long n_element = problem_pt->mesh_pt()->nelement();
  for(unsigned long e = 0;e<n_element;e++)
   {
    GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the ndofs in each element
    unsigned n_var = elem_pt->ndof();
    //Get the jacobian matrix
    DenseMatrix<double> jac(n_var);
    // storage for the residual
    Vector<double> res(n_var);
    //Get unperturbed raw jacobian
    elem_pt->get_jacobian(res,jac);

    //Multiply the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      for(unsigned m=0;m<n_var;m++)
       {
        unsigned unknown = elem_pt->eqn_number(m);
        Jy[eqn_number] += jac(n,m)*y[unknown];
       }
     }
   }
  //The final entry of the vector will be zero

  //Now resolve to find alpha
  Linear_solver_pt->resolve(Jy,*Alpha_pt);

  //The vector that multiplie the product matrix is actually y - alpha
  DoubleVector y_minus_alpha(this->distribution_pt(),0.0);
  for(unsigned n=0;n<n_dof;n++)
   {
    y_minus_alpha[n] = y[n] - (*Alpha_pt)[n];
   }

  //We can now construct our multipliers
  //Prepare to scale
  double dof_length=0.0, a_length=0.0, alpha_length=0.0;
  for(unsigned n=0;n<n_dof;n++)
   {
    if(std::fabs(problem_pt->dof(n)) > dof_length)
     {dof_length = std::fabs(problem_pt->dof(n));}
    if(std::fabs(a[n]) > a_length) {a_length = std::fabs(a[n]);}
    if(std::fabs(y_minus_alpha[n]) > alpha_length)
     {alpha_length = std::fabs(y_minus_alpha[n]);}
   }

  double a_mult = dof_length/a_length;
  double alpha_mult = dof_length/alpha_length;
  const double FD_step = 1.0e-8;
  a_mult += FD_step; alpha_mult += FD_step;
  a_mult *= FD_step; alpha_mult *= FD_step;

  //Local storage for the product terms
  DoubleVector Jprod_a(this->distribution_pt(),0.0),
   Jprod_alpha(this->distribution_pt(),0.0);

  //Calculate the product of the jacobian matrices, etc
  for(unsigned long e = 0;e<n_element;e++)
   {
    GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the ndofs in each element
    unsigned n_var = handler_pt->ndof(elem_pt);
    //Get the jacobian matrices
    DenseMatrix<double> jac(n_var), jac_a(n_var), jac_alpha(n_var);
    // elemental residual storage
    Vector<double> res_elemental(n_var);
    //Get unperturbed jacobian
    handler_pt->get_jacobian(elem_pt,res_elemental,jac);

    //Backup the dofs
    Vector<double> dof_bac(n_var);
    //Perturb the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
      dof_bac[n] = problem_pt->dof(eqn_number);
      //Pertub by vector a
      problem_pt->dof(eqn_number) += a_mult*a[eqn_number];
     }

    problem_pt->actions_after_change_in_bifurcation_parameter();

    //Now get the new jacobian
    handler_pt->get_jacobian(elem_pt,res_elemental,jac_a);

    //Perturb the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
      problem_pt->dof(eqn_number) = dof_bac[n];
      //Pertub by vector a
      problem_pt->dof(eqn_number) += alpha_mult*y_minus_alpha[eqn_number];
     }

    problem_pt->actions_after_change_in_bifurcation_parameter();

    //Now get the new jacobian
    handler_pt->get_jacobian(elem_pt,res_elemental,jac_alpha);

    //Reset the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
      problem_pt->dof(eqn_number) = dof_bac[n];
     }

    problem_pt->actions_after_change_in_bifurcation_parameter();

    //OK, now work out the products
    //Note the (n_var-1), we are only interested in the non-augmented
    //jacobian
    for(unsigned n=0;n<(n_var-1);n++)
     {
      unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
      double prod_a=0.0, prod_alpha=0.0;
      for(unsigned m=0;m<(n_var-1);m++)
       {
        unsigned unknown = handler_pt->eqn_number(elem_pt,m);
        prod_a += (jac_a(n,m) - jac(n,m))*y[unknown];
        prod_alpha += (jac_alpha(n,m) - jac(n,m))*y[unknown];
       }
      Jprod_a[eqn_number] += prod_a/a_mult;
      Jprod_alpha[eqn_number] += prod_alpha/alpha_mult;
     }
   }

   Jprod_alpha[n_dof-1] = 0.0;
   Jprod_a[n_dof-1] = 0.0;

   //OK, now we can formulate the next vectors
   //The assumption here is that the result has been set to the
   //residuals.
   for(unsigned n=0;n<n_dof-1;n++)
    {
     b[n] = result[n_dof+n] - Jprod_a[n];
    }
   //The final residual is the entry corresponding to the original parameter
   b[n_dof-1] = result[n_dof-1];

   //Allocate storage for E which can be used in the resolve
   if(E_pt!=0) {delete E_pt;}
   E_pt = new DoubleVector(this->distribution_pt(),0.0);
   DoubleVector f(this->distribution_pt(),0.0);
   Linear_solver_pt->resolve(b,f);
   Linear_solver_pt->resolve(Jprod_alpha,*E_pt);

   //Calculate the final entry in the vector e
   const double e_final = (*E_pt)[n_dof-1];
   //Calculate the final entry in the vector d
   const double d_final = f[n_dof-1]/e_final;
   //Assemble the final corrections
   for(unsigned n=0;n<n_dof-1;n++)
    {
     result[n] = a[n] - (*Alpha_pt)[n]*d_final + d_final*y[n];
     result[n_dof+n] = f[n] - (*E_pt)[n]*d_final;
    }
   //The result corresponding to the parameter
   result[n_dof-1] = a[n_dof-1] - (*Alpha_pt)[n_dof-1]*d_final;

   //The sign of the jacobian is the sign of the final entry in e
   problem_pt->sign_of_jacobian() =
    static_cast<int>(std::fabs(e_final)/e_final);

   //Switch things to our block solver
   handler_pt->solve_full_system();

   //If we are not storing things, clear the results
   if(!Enable_resolve)
    {
     //We no longer need to store the matrix
     Linear_solver_pt->disable_resolve();
     delete Alpha_pt; Alpha_pt=0;
     delete E_pt; E_pt=0;
    }
   //Otherwise, also store the problem pointer
   else
    {
     Problem_pt = problem_pt;
    }
 }


 //======================================================================
 //Hack the re-solve to use the block factorisation
 //======================================================================
 void AugmentedBlockFoldLinearSolver::resolve(const DoubleVector &rhs,
                                              DoubleVector &result)
 {
  //Check that the factors have been stored
  if(Alpha_pt==0)
   {
    throw OomphLibError("The required vectors have not been stored",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

  //Get the pointer to the problem
  Problem* const problem_pt = Problem_pt;

  FoldHandler* handler_pt =
   static_cast<FoldHandler*>(problem_pt->assembly_handler_pt());

  //Switch things to our block solver
  handler_pt->solve_augmented_block_system();

  //We need to find out the number of dofs in the problem
  unsigned n_dof = problem_pt->ndof();

#ifdef PARANOID
  // if the result is setup then it should not be distributed
  if (result.built())
   {
    if (result.distributed())
     {
      throw OomphLibError("The result vector must not be distributed",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
  // the rhs must be setup
  if (!rhs.built())
   {
      throw OomphLibError("The rhs vector must be setup",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
   }
#endif

  // if the result vector is not setup then rebuild with distribution = global
  if (!result.built())
   {
    result.build(rhs.distribution_pt(),0.0);
   }

  //Setup storage
  DoubleVector a(this->distribution_pt(),0.0), b(this->distribution_pt(),0.0);

  //Set the values of the a vector
  for(unsigned n=0;n<(n_dof-1);n++) {a[n] = rhs[n];}
  //The entry associated with the additional parameter is zero
  a[n_dof-1] = 0.0;

  Linear_solver_pt->enable_resolve();

  // Copy rhs vector into local storage so it doesn't get overwritten
  // if the linear solver decides to initialise the solution vector, say,
  // which it's quite entitled to do!
  DoubleVector input_a(a);

  Linear_solver_pt->resolve(input_a,a);

  //We can now construct our multipliers
  //Prepare to scale
  double dof_length=0.0, a_length=0.0;
  for(unsigned n=0;n<n_dof;n++)
   {
    if(std::fabs(problem_pt->dof(n)) > dof_length)
     {dof_length = std::fabs(problem_pt->dof(n));}

    if(std::fabs(a[n]) > a_length) {a_length = std::fabs(a[n]);}
   }
  double a_mult = dof_length/a_length;
  const double FD_step = 1.0e-8;
  a_mult += FD_step;
  a_mult *= FD_step;

  DoubleVector Jprod_a(this->distribution_pt(),0.0);

  unsigned long n_element = problem_pt->mesh_pt()->nelement();
   for(unsigned long e = 0;e<n_element;e++)
    {
     GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
     //Loop over the ndofs in each element
     unsigned n_var = handler_pt->ndof(elem_pt);
     //Get some jacobian matrices
     DenseMatrix<double> jac(n_var), jac_a(n_var);
     // the elemental residual
     Vector<double> res_elemental(n_var);
     //Get unperturbed jacobian
     handler_pt->get_jacobian(elem_pt,res_elemental,jac);

     //Backup the dofs
     Vector<double> dof_bac(n_var);
     //Perturb the dofs
     for(unsigned n=0;n<n_var;n++)
      {
       unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
       dof_bac[n] = problem_pt->dof(eqn_number);
       //Pertub by vector a
       problem_pt->dof(eqn_number) += a_mult*a[eqn_number];
      }

     problem_pt->actions_after_change_in_bifurcation_parameter();

     //Now get the new jacobian
     handler_pt->get_jacobian(elem_pt,res_elemental,jac_a);

     //Reset the dofs
     for(unsigned n=0;n<n_var;n++)
      {
       unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
       problem_pt->dof(eqn_number) = dof_bac[n];
      }

     problem_pt->actions_after_change_in_bifurcation_parameter();

     //OK, now work out the products
     for(unsigned n=0;n<(n_var-1);n++)
      {
       unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
       double prod_a=0.0;
       for(unsigned m=0;m<(n_var-1);m++)
        {
         unsigned unknown = handler_pt->eqn_number(elem_pt,m);
         prod_a += (jac_a(n,m) - jac(n,m))*handler_pt->Y[unknown];
        }
       Jprod_a[eqn_number] += prod_a/a_mult;
      }
    }

   Jprod_a[n_dof-1] = 0.0;

   //OK, now we can formulate the next vectors
   for(unsigned n=0;n<n_dof-1;n++)
    {
     b[n] = rhs[n_dof+n] - Jprod_a[n];
    }
   //The residuals for the final entry should be those associated
   //with the parameter
   b[n_dof-1] = rhs[n_dof-1];

   DoubleVector f(this->distribution_pt(),0.0);

   Linear_solver_pt->resolve(b,f);

   //Calculate the final entry in the vector d
   const double d_final = f[n_dof-1]/(*E_pt)[n_dof-1];
   //Assemble the final corrections
   for(unsigned n=0;n<n_dof-1;n++)
    {
     result[n] = a[n] - (*Alpha_pt)[n]*d_final + d_final*handler_pt->Y[n];
     result[n_dof+n] = f[n] - (*E_pt)[n]*d_final;
    }
   //The result corresponding to the paramater
   result[n_dof-1] = a[n_dof-1] - (*Alpha_pt)[n_dof-1]*d_final;

   Linear_solver_pt->disable_resolve();

   //Switch things to our block solver
   handler_pt->solve_full_system();
  }


 ///////////////////////////////////////////////////////////////////////
 // Non-inline functions for the FoldHandler class
 //////////////////////////////////////////////////////////////////////


 //========================================================================
 /// Constructor: Initialise the fold handler,
 /// by setting initial guesses for Y, Phi and calculating Count.
 /// If the system changes, a new  handler must be constructed.
 //========================================================================
 FoldHandler::FoldHandler(Problem* const &problem_pt,
                          double* const &parameter_pt) :
  Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
 {
  //Set the problem pointer
  Problem_pt = problem_pt;
  //Set the number of degrees of freedom
  Ndof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  LinearAlgebraDistribution* dist_pt = new
   LinearAlgebraDistribution(problem_pt->communicator_pt(),Ndof,false);

  //Resize the vectors of additional dofs and constants
  Phi.resize(Ndof);
  Y.resize(Ndof);
  Count.resize(Ndof,0);

  //Loop over all the elements in the problem
  unsigned n_element = problem_pt->mesh_pt()->nelement();
  for(unsigned e=0;e<n_element;e++)
   {
    GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the local freedoms in the element
    unsigned n_var = elem_pt->ndof();
    for(unsigned n=0;n<n_var;n++)
     {
      //Increase the associated global equation number counter
      ++Count[elem_pt->eqn_number(n)];
     }
   }

  //Calculate the value Phi by
  //solving the system JPhi = dF/dlambda

  //Locally cache the linear solver
  LinearSolver* const linear_solver_pt = problem_pt->linear_solver_pt();

  //Save the status before entry to this routine
  bool enable_resolve = linear_solver_pt->is_resolve_enabled();

  //We need to do a resolve
  linear_solver_pt->enable_resolve();

  //Storage for the solution
  DoubleVector x(dist_pt,0.0);

  //Solve the standard problem, we only want to make sure that
  //we factorise the matrix, if it has not been factorised. We shall
  //ignore the return value of x.
  linear_solver_pt->solve(problem_pt,x);

  //Get the vector dresiduals/dparameter
  problem_pt->get_derivative_wrt_global_parameter(parameter_pt,x);

  // Copy rhs vector into local storage so it doesn't get overwritten
  // if the linear solver decides to initialise the solution vector, say,
  // which it's quite entitled to do!
  DoubleVector input_x(x);

  //Now resolve the system with the new RHS and overwrite the solution
  linear_solver_pt->resolve(input_x,x);

  //Restore the storage status of the linear solver
  if (enable_resolve)
   {
    linear_solver_pt->enable_resolve();
   }
  else
   {
    linear_solver_pt->disable_resolve();
   }

  //Add the global parameter as an unknown to the problem
  problem_pt->Dof_pt.push_back(parameter_pt);


  //Normalise the initial guesses for phi
  double length = 0.0;
  for(unsigned n=0;n<Ndof;n++) {length += x[n]*x[n];}
  length = sqrt(length);

  //Now add the null space components to the problem unknowns
  //and initialise them and Phi to the same normalised values
  for(unsigned n=0;n<Ndof;n++)
   {
    problem_pt->Dof_pt.push_back(&Y[n]);
    Y[n] = Phi[n] = -x[n]/length;
   }

  // delete the dist_pt
  problem_pt->Dof_distribution_pt->build(problem_pt->communicator_pt(),
                                         Ndof*2+1,true);
  //Remove all previous sparse storage used during Jacobian assembly
  Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);

  delete dist_pt;
 }


 ///Constructor in which the eigenvector is passed as an initial
 ///guess
 FoldHandler::FoldHandler(Problem* const &problem_pt,
                          double* const &parameter_pt,
                          const DoubleVector &eigenvector) :
  Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
 {
  //Set the problem pointer
  Problem_pt = problem_pt;
  //Set the number of degrees of freedom
  Ndof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  LinearAlgebraDistribution* dist_pt = new
   LinearAlgebraDistribution(problem_pt->communicator_pt(),
                             Ndof,false);

  //Resize the vectors of additional dofs and constants
  Phi.resize(Ndof);
  Y.resize(Ndof);
  Count.resize(Ndof,0);

  //Loop over all the elements in the problem
  unsigned n_element = problem_pt->mesh_pt()->nelement();
  for(unsigned e=0;e<n_element;e++)
   {
    GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the local freedoms in the element
    unsigned n_var = elem_pt->ndof();
    for(unsigned n=0;n<n_var;n++)
     {
      //Increase the associated global equation number counter
      ++Count[elem_pt->eqn_number(n)];
     }
   }

  //Add the global parameter as an unknown to the problem
  problem_pt->Dof_pt.push_back(parameter_pt);

  
  //Normalise the initial guesses for the eigenvecor
  double length = 0.0;
  for(unsigned n=0;n<Ndof;n++)
   {length += eigenvector[n]*eigenvector[n];}
  length = sqrt(length);

  //Now add the null space components to the problem unknowns
  //and initialise them and Phi to the same normalised values
  for(unsigned n=0;n<Ndof;n++)
   {
    problem_pt->Dof_pt.push_back(&Y[n]);
    Y[n] = Phi[n] = eigenvector[n]/length;
   }

  // delete the dist_pt
  problem_pt->Dof_distribution_pt->build(problem_pt->communicator_pt(),
                                         Ndof*2+1,true);
  //Remove all previous sparse storage used during Jacobian assembly
  Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);

  delete dist_pt;
 }


  ///Constructor in which the eigenvector and normalisation
  ///vector are  passed as an initial guess
 FoldHandler::FoldHandler(Problem* const &problem_pt,
                          double* const &parameter_pt,
                          const DoubleVector &eigenvector,
                          const DoubleVector &normalisation) :
  Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
 {
  //Set the problem pointer
  Problem_pt = problem_pt;
  //Set the number of degrees of freedom
  Ndof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  LinearAlgebraDistribution* dist_pt = new
   LinearAlgebraDistribution(problem_pt->communicator_pt(),
                             Ndof,false);

  //Resize the vectors of additional dofs and constants
  Phi.resize(Ndof);
  Y.resize(Ndof);
  Count.resize(Ndof,0);

  //Loop over all the elements in the problem
  unsigned n_element = problem_pt->mesh_pt()->nelement();
  for(unsigned e=0;e<n_element;e++)
   {
    GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the local freedoms in the element
    unsigned n_var = elem_pt->ndof();
    for(unsigned n=0;n<n_var;n++)
     {
      //Increase the associated global equation number counter
      ++Count[elem_pt->eqn_number(n)];
     }
   }

  //Add the global parameter as an unknown to the problem
  problem_pt->Dof_pt.push_back(parameter_pt);

  
  //Normalise the initial guesses for the eigenvecor
  double length = 0.0;
  for(unsigned n=0;n<Ndof;n++)
   {length += eigenvector[n]*normalisation[n];}
  length = sqrt(length);

  //Now add the null space components to the problem unknowns
  //and initialise them and Phi to the same normalised values
  for(unsigned n=0;n<Ndof;n++)
   {
    problem_pt->Dof_pt.push_back(&Y[n]);
    Y[n] = eigenvector[n]/length;
    Phi[n] = normalisation[n];
   }

  // delete the dist_pt
  problem_pt->Dof_distribution_pt->build(problem_pt->communicator_pt(),
                                         Ndof*2+1,true);
  //Remove all previous sparse storage used during Jacobian assembly
  Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);

  delete dist_pt;
 }


 
 //=================================================================
 /// The number of unknowns is 2n+1
 //================================================================
 unsigned FoldHandler::ndof(GeneralisedElement* const &elem_pt)
 {
  unsigned raw_ndof = elem_pt->ndof();
  //Return different values depending on the type of block decomposition
  switch(Solve_which_system)
   {
   case Full_augmented:
    return (2*raw_ndof + 1);
    break;

   case Block_augmented_J:
    return (raw_ndof+1);
    break;

   case Block_J:
    return raw_ndof;
    break;

   default:
    std::ostringstream error_stream;
    error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
                 << " not " << Solve_which_system << "\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 }

 //=====================================================================
 /// Return the global equation number associated with local equation
 /// number ieqn_local. We choose to number the unknowns according
 /// to the augmented system.
 //=======================================================================
 unsigned long FoldHandler::eqn_number(GeneralisedElement* const &elem_pt,
                                       const unsigned &ieqn_local)
 {
  //Find the number of non-augmented dofs in the element
  unsigned raw_ndof = elem_pt->ndof();
  //Storage for the global eqn number
  unsigned long global_eqn=0;
  //If we are a "standard" unknown, just return the standard equation number
  if(ieqn_local < raw_ndof)
   {
    global_eqn = elem_pt->eqn_number(ieqn_local);
   }
  //Otherwise if we are at an unknown corresponding to the bifurcation
  //parameter return the global equation number of the parameter
  else if(ieqn_local==raw_ndof)
   {
    global_eqn = Ndof;
   }
  //Otherwise we are in the unknown corresponding to a null vector
  //return the global unknown Ndof + 1 + global unknown of "standard" unknown.
  else
   {
    global_eqn = Ndof + 1 + elem_pt->eqn_number(ieqn_local - 1 - raw_ndof);
   }

  //Return the global equation number
  return global_eqn;
 }

 //====================================================================
 ///Formulate the augmented system
 //====================================================================
 void FoldHandler::get_residuals(GeneralisedElement* const &elem_pt,
                                 Vector<double> &residuals)
 {
  //Need to get raw residuals and jacobian
  unsigned raw_ndof = elem_pt->ndof();

  //Find out which system we are solving
  switch(Solve_which_system)
   {
    //If we are solving the standard system
   case Block_J:
   {
    //Get the basic residuals
    elem_pt->get_residuals(residuals);
   }
   break;

   //If we are solving the augmented-by-one system
   case Block_augmented_J:
   {
    //Get the basic residuals
    elem_pt->get_residuals(residuals);

    //Zero the final residual
    residuals[raw_ndof] = 0.0;
   }
   break;

   //If we are solving the full augmented system
   case Full_augmented:
   {
    DenseMatrix<double> jacobian(raw_ndof);
    //Get the basic residuals and jacobian initially
    elem_pt->get_jacobian(residuals,jacobian);

    //The normalisation equation must be initialised to -1.0/number of elements
    residuals[raw_ndof] = -1.0/Problem_pt->mesh_pt()->nelement();

    //Now assemble the equations Jy = 0
    for(unsigned i=0;i<raw_ndof;i++)
     {
      residuals[raw_ndof+1+i] = 0.0;
      for(unsigned j=0;j<raw_ndof;j++)
       {
        residuals[raw_ndof+1+i] += jacobian(i,j)*Y[elem_pt->eqn_number(j)];
       }
      //Add the contribution to phi.y=1
      //Need to divide by the number of elements that contribute to this
      //unknown so that we do get phi.y exactly.
      unsigned global_eqn = elem_pt->eqn_number(i);
      residuals[raw_ndof] +=
       (Phi[global_eqn]*Y[global_eqn])/Count[global_eqn];
     }
   }
   break;

   default:
    std::ostringstream error_stream;
    error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
                 << " not " << Solve_which_system << "\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 }


 //=================================================================
 /// Calculate the elemental Jacobian matrix "d equation
 /// / d variable" for the augmented system
 //=================================================================
 void FoldHandler::get_jacobian(GeneralisedElement* const &elem_pt,
                                Vector<double> &residuals,
                                DenseMatrix<double> &jacobian)
 {

  //Find the number of augmented dofs
  unsigned augmented_ndof = ndof(elem_pt);
  //Find the non-augmented dofs
  unsigned raw_ndof = elem_pt->ndof();

  //Which system are we solving
  switch(Solve_which_system)
   {
    //If we are solving the original system
   case Block_J:
   {
    //Just get the raw jacobian and residuals
    elem_pt->get_jacobian(residuals,jacobian);
   }
   break;

   //If we are solving the augmented-by-one system
   case Block_augmented_J:
   {
    //Get the full residuals, we need them
    get_residuals(elem_pt,residuals);

    //Need to get the raw jacobian (and raw residuals)
    Vector<double> newres(augmented_ndof);
    elem_pt->get_jacobian(newres,jacobian);

    //Now do finite differencing stuff
    const double FD_step = 1.0e-8;
    //Fill in the first lot of finite differences
    {
     //increase the global parameter
     double *unknown_pt = Problem_pt->Dof_pt[Ndof];
     double init = *unknown_pt;
     *unknown_pt += FD_step;

     //Now do any possible updates
     Problem_pt
      ->actions_after_change_in_bifurcation_parameter();

     //Get the new (modified) residuals
     get_residuals(elem_pt,newres);

     //The final column  is given by the difference
     //between the residuals
     for(unsigned n=0;n<raw_ndof;n++)
      {
       jacobian(n,augmented_ndof-1) = (newres[n] - residuals[n])/FD_step;
      }
     //Reset the global parameter
     *unknown_pt = init;

     //Now do any possible updates
     Problem_pt
      ->actions_after_change_in_bifurcation_parameter();
    }

    //Fill in the bottom row
    for(unsigned n=0;n<raw_ndof;n++)
     {
      unsigned local_eqn = elem_pt->eqn_number(n);
      jacobian(augmented_ndof-1,n)
       = Phi[local_eqn]/Count[local_eqn];
     }
   }
   break;

  //Otherwise solving the full system
   case Full_augmented:
   {
    //Get the residuals for the augmented system
    get_residuals(elem_pt,residuals);

    //Need to get the raw residuals and jacobian
    Vector<double> newres(raw_ndof);
    DenseMatrix<double> newjac(raw_ndof);
    elem_pt->get_jacobian(newres,jacobian);

    //Fill in the jacobian on the diagonal sub-block of
    //the null-space equations
    for(unsigned n=0;n<raw_ndof;n++)
     {
      for(unsigned m=0;m<raw_ndof;m++)
       {
        jacobian(raw_ndof+1+n,raw_ndof+1+m) = jacobian(n,m);
       }
     }

    //Now finite difference wrt the global unknown
    const double FD_step = 1.0e-8;
    //Fill in the first lot of finite differences
    {
   //increase the global parameter
     double *unknown_pt = Problem_pt->Dof_pt[Ndof];
     double init = *unknown_pt;
     *unknown_pt += FD_step;
     //Need to update the function
     Problem_pt->actions_after_change_in_bifurcation_parameter();

     //Get the new raw residuals and jacobian
     elem_pt->get_jacobian(newres,newjac);

     //The end of the first row is given by the difference
     //between the residuals
     for(unsigned n=0;n<raw_ndof;n++)
      {
       jacobian(n,raw_ndof) = (newres[n] - residuals[n])/FD_step;
       //The end of the second row is given by the difference multiplied by
       //the product
       for(unsigned l=0;l<raw_ndof;l++)
        {
         jacobian(raw_ndof+1+n,raw_ndof)
          += (newjac(n,l) - jacobian(n,l))*Y[elem_pt->eqn_number(l)]/
          FD_step;
        }
      }
     //Reset the global parameter
     *unknown_pt = init;

     //Need to update the function
     Problem_pt->actions_after_change_in_bifurcation_parameter();
    }

    //Now fill in the first column of the second rows
    {
     for(unsigned n=0;n<raw_ndof;n++)
      {
       unsigned long global_eqn = eqn_number(elem_pt,n);
       //Increase the first lot
       double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
       double init = *unknown_pt;
       *unknown_pt += FD_step;
       Problem_pt->actions_before_newton_convergence_check();/// ALICE

       //Get the new jacobian
       elem_pt->get_jacobian(newres,newjac);

       //Work out the differences
       for(unsigned k=0;k<raw_ndof;k++)
        {
         //jacobian(raw_ndof+k,n) = 0.0;
         for(unsigned l=0;l<raw_ndof;l++)
          {
           jacobian(raw_ndof+1+k,n)
            += (newjac(k,l) - jacobian(k,l))*Y[elem_pt->eqn_number(l)]/
            FD_step;
          }
        }
       *unknown_pt = init;
       Problem_pt->actions_before_newton_convergence_check();/// ALICE
      }
    }

    //Fill in the row corresponding to the parameter
    for(unsigned n=0;n<raw_ndof;n++)
     {
      unsigned global_eqn = elem_pt->eqn_number(n);
      jacobian(raw_ndof,raw_ndof+1+n)
       = Phi[global_eqn]/Count[global_eqn];
     }

  //   //Loop over the dofs
//     for(unsigned n=0;n<augmented_ndof;n++)
//      {
//       unsigned long global_eqn = eqn_number(elem_pt,n);
//       double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
//       double init = *unknown_pt;
//       *unknown_pt += FD_step;

//       //Get the new residuals
//       get_residuals(elem_pt,newres);

//       for(unsigned m=0;m<augmented_ndof;m++)
//        {
//         jacobian(m,n) = (newres[m] - residuals[m])/FD_step;
//        }
//       //Reset the unknown
//       *unknown_pt = init;
//      }
   }
   break;

   default:
    std::ostringstream error_stream;
    error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
                 << " not " << Solve_which_system << "\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);

   }
 }



 //====================================================================
 ///Formulate the derivatives of the augmented system with respect
 ///to a parameter
 //====================================================================
 void FoldHandler::get_dresiduals_dparameter(
  GeneralisedElement* const &elem_pt,
  double* const &parameter_pt,
  Vector<double> &dres_dparam)
 {
  //Need to get raw residuals and jacobian
  unsigned raw_ndof = elem_pt->ndof();

  //Find out which system we are solving
  switch(Solve_which_system)
   {
    //If we are solving the standard system
   case Block_J:
   {
    //Get the basic residual derivatives
    elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);
   }
   break;

   //If we are solving the augmented-by-one system
   case Block_augmented_J:
   {
    //Get the basic residual derivatives
    elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);

    //Zero the final derivative
    dres_dparam[raw_ndof] = 0.0;
   }
   break;

   //If we are solving the full augmented system
   case Full_augmented:
   {
    DenseMatrix<double> djac_dparam(raw_ndof);
    //Get the basic residuals and jacobian derivatives initially
    elem_pt->get_djacobian_dparameter(parameter_pt,dres_dparam,djac_dparam);

    //The normalisation equation does not depend on the parameter
    dres_dparam[raw_ndof] = 0.0;

    //Now assemble the equations dJy/dparameter = 0
    for(unsigned i=0;i<raw_ndof;i++)
     {
      dres_dparam[raw_ndof+1+i] = 0.0;
      for(unsigned j=0;j<raw_ndof;j++)
       {
        dres_dparam[raw_ndof+1+i] += djac_dparam(i,j)*
         Y[elem_pt->eqn_number(j)];
       }
     }
   }
   break;

   default:
    std::ostringstream error_stream;
    error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
                 << " not " << Solve_which_system << "\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 }


 //========================================================================
 /// Overload the derivative of the residuals and jacobian
 /// with respect to a parameter so that it breaks because it should not
 /// be required
 //========================================================================
 void FoldHandler::get_djacobian_dparameter(GeneralisedElement* const &elem_pt,
                                            double* const &parameter_pt,
                                            Vector<double> &dres_dparam,
                                            DenseMatrix<double> &djac_dparam)
 {
  std::ostringstream error_stream;
  error_stream <<
   "This function has not been implemented because it is not required\n";
  error_stream << "in standard problems.\n";
  error_stream <<
   "If you find that you need it, you will have to implement it!\n\n";

  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
 }


 //=====================================================================
 /// Overload the hessian vector product function so that
 /// it breaks because it should not be required
 //========================================================================
 void FoldHandler::get_hessian_vector_products(
  GeneralisedElement* const &elem_pt,
  Vector<double> const &Y,
  DenseMatrix<double> const &C,
  DenseMatrix<double> &product)
 {
  std::ostringstream error_stream;
  error_stream <<
   "This function has not been implemented because it is not required\n";
  error_stream << "in standard problems.\n";
  error_stream <<
   "If you find that you need it, you will have to implement it!\n\n";

  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
 }


 //==========================================================================
 /// Return the eigenfunction(s) associated with the bifurcation that
 /// has been detected in bifurcation tracking problems
 //==========================================================================
 void FoldHandler::get_eigenfunction(Vector<DoubleVector> &eigenfunction)
 {
  //There is only one (real) null vector
  eigenfunction.resize(1);
  LinearAlgebraDistribution dist(Problem_pt->communicator_pt(),Ndof,false);
  //Rebuild the vector
  eigenfunction[0].build(&dist,0.0);
  //Set the value to be the null vector
  for(unsigned n=0;n<Ndof;n++)
   {
    eigenfunction[0][n] = Y[n];
   }
 }



 //=======================================================================
 /// Destructor return the problem to its original (non-augmented) state
 //=======================================================================
 FoldHandler::~FoldHandler()
 {
  //If we are using the block solver reset the problem's linear solver
  //to the original one
  AugmentedBlockFoldLinearSolver *block_fold_solver_pt
   = dynamic_cast<AugmentedBlockFoldLinearSolver*>(
    Problem_pt->linear_solver_pt());

  if(block_fold_solver_pt)
   {
    //Reset the problem's linear solver
    Problem_pt->linear_solver_pt() = block_fold_solver_pt->linear_solver_pt();
    //Delete the block solver
    delete block_fold_solver_pt;
   }

  //Resize the number of dofs
  Problem_pt->Dof_pt.resize(Ndof);
  Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                         Ndof,false);
  //Remove all previous sparse storage used during Jacobian assembly
  Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 }

 //====================================================================
 /// Set to solve the augmented-by-one block system.
 //===================================================================
 void FoldHandler::solve_augmented_block_system()
 {
  //Only bother to do anything if we haven't already set the flag
  if(Solve_which_system!=Block_augmented_J)
   {
    //If we were solving the system with the original jacobian add the
    //parameter
    if(Solve_which_system==Block_J)
     {
      Problem_pt->Dof_pt.push_back(Parameter_pt);
     }

    //Restrict the problem to the standard variables and
    //the bifurcation parameter only
    Problem_pt->Dof_pt.resize(Ndof+1);
    Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                           Ndof+1,false);
    //Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
    //Set the solve flag
    Solve_which_system = Block_augmented_J;
   }
 }


 //====================================================================
 /// Set to solve the block system. The boolean flag specifies
 /// whether the block decomposition should use exactly the same jacobian
 //===================================================================
 void FoldHandler::solve_block_system()
 {
  //Only bother to do anything if we haven't already set the flag
  if(Solve_which_system!=Block_J)
   {
    //Restrict the problem to the standard variables
    Problem_pt->Dof_pt.resize(Ndof);
    Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                           Ndof,false);
    //Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);

    //Set the solve flag
    Solve_which_system = Block_J;
   }
 }

 //=================================================================
 /// Set to Solve non-block system
 //=================================================================
 void FoldHandler::solve_full_system()
 {
  //Only do something if we are not solving the full system
  if(Solve_which_system!=Full_augmented)
   {
    //If we were solving the problem without any augmentation,
    //add the parameter again
    if(Solve_which_system==Block_J)
     {Problem_pt->Dof_pt.push_back(Parameter_pt);}

    //Always add the additional unknowns back into the problem
    for(unsigned n=0;n<Ndof;n++)
     {
      Problem_pt->Dof_pt.push_back(&Y[n]);
     }

    // update the Dof distribution pt
    Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                           Ndof*2+1,false);
    //Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);

    Solve_which_system = Full_augmented;
   }
 }



 //======================================================================
 /// Clean up the memory that may have been allocated by the solver
 //=====================================================================
 BlockPitchForkLinearSolver::~BlockPitchForkLinearSolver()
 {
  if(B_pt!=0) {delete B_pt;}
  if(C_pt!=0) {delete C_pt;}
  if(D_pt!=0) {delete D_pt;}
  if(dJy_dparam_pt!=0) {delete dJy_dparam_pt;}
 }

 //===================================================================
 /// Use a block factorisation to solve the augmented system
 /// associated with a PitchFork bifurcation.
 //===================================================================
 void BlockPitchForkLinearSolver:: solve(Problem* const &problem_pt,
                                         DoubleVector &result)
 {
     std::cout << "Block pitchfork solve" << std::endl;
  //Locally cache the pointer to the handler.
  PitchForkHandler* handler_pt =
   static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());

  //Get the augmented distribution from the handler and use it as
  //the distribution of the linear solver
  LinearAlgebraDistribution aug_dist =
   handler_pt->Augmented_dof_distribution_pt;
  this->build_distribution(aug_dist);

  // If the result vector is not setup then die
  if(!result.built())
   {
    throw OomphLibError("Result vector must be built\n",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

  //Store the distribution of the result, which is probably not in
  //the natural distribution of the augmented solver
  LinearAlgebraDistribution result_dist(result.distribution_pt());

  //Locally cache a pointer to the parameter
  double* const parameter_pt = handler_pt->Parameter_pt;

  //Firstly get the derivatives of the full augmented system with
  //respect to the parameter

  //Allocate storage for the derivatives of the residuals with respect
  //to the global parameter using the augmented distribution
  DoubleVector dRdparam;
  //Then get the appropriate derivatives which will come back in
  //some distribution
  problem_pt->get_derivative_wrt_global_parameter(parameter_pt,dRdparam);
  //Redistribute into the augmented distribution
  dRdparam.redistribute(&aug_dist);

  //Switch the handler to "block solver" mode (sort out distribution)
  handler_pt->solve_block_system();

  //Temporary vector for the result (I shouldn't have to set this up)
  DoubleVector x1;

  //We are going to do resolves using the underlying linear solver
  Linear_solver_pt->enable_resolve();
  //Solve the first (standard) system Jx1 = R
  Linear_solver_pt->solve(problem_pt,x1);

  //Allocate storage for B, C and D which can be used in the resolve
  //and the derivatives of the jacobian/eigenvector product with
  //respect to the parameter
  if(B_pt!=0) {delete B_pt;}
  B_pt = new DoubleVector(Linear_solver_pt->distribution_pt(),0.0);
  if(C_pt!=0) {delete C_pt;}
  C_pt = new DoubleVector(Linear_solver_pt->distribution_pt(),0.0);
  if(D_pt!=0) {delete D_pt;}
  D_pt = new DoubleVector(Linear_solver_pt->distribution_pt(),0.0);
  //Need this to be in the distribution of the dofs
  if(dJy_dparam_pt!=0) {delete dJy_dparam_pt;}
  dJy_dparam_pt = new DoubleVector(handler_pt->Dof_distribution_pt,0.0);

  //Get the symmetry vector from the handler should have the
  //Dof distribution
  DoubleVector psi = handler_pt->Psi;

  //Temporary vector for the rhs that is dR/dparam (augmented distribution)
  //DoubleVector F(Linear_solver_pt->distribution_pt(),0.0);
  DoubleVector F(handler_pt->Dof_distribution_pt);
  const unsigned n_dof_local = F.nrow_local();

  //f.nrow local copied from dRdparam
  for(unsigned n=0;n<n_dof_local;n++) {F[n] = dRdparam[n];}
  //Fill in the rhs that is dJy/dparam  //dJy_dparam nrow local

  //In standard cases the offset here will be one
  unsigned offset=1;
#ifdef OOMPH_HAS_MPI
  //If we are distributed and not on the first processor
  //then there is no offset and the eigenfunction is
  //directly after the standard dofs
  int my_rank = problem_pt->communicator_pt()->my_rank();
  if((problem_pt->distributed()) && (my_rank!=0)) {offset=0;}
#endif
  for(unsigned n=0;n<n_dof_local;n++)
   {(*dJy_dparam_pt)[n] = dRdparam[n_dof_local + offset + n];}

  //Now resolve to find c and d
  //First, redistribute F and psi into the linear algebra distribution
  F.redistribute(Linear_solver_pt->distribution_pt());
  psi.redistribute(Linear_solver_pt->distribution_pt());

  Linear_solver_pt->resolve(F,*C_pt);
  Linear_solver_pt->resolve(psi,*D_pt);

  //We can now construct various dot products
  double psi_d = psi.dot(*D_pt);
  double psi_c = psi.dot(*C_pt);
  double psi_x1 = psi.dot(x1);

  //Calculate another intermediate constant
  const double Psi = psi_d/psi_c;

  //Construct the second intermediate value,
  //assumption is that result has been
  //set to the current value of the residuals
  //First, redistribute into the Natural distribution of the augmented system
  result.redistribute(&aug_dist);
  //The required parameter is that corresponding to the dof, which
  //is only stored on the root processor
  double parameter_residual = result[n_dof_local];
#ifdef OOMPH_HAS_MPI
  //Broadcast to all others, if we have a distributed problem
  if(problem_pt->distributed())
   {
    MPI_Bcast(&parameter_residual,1,MPI_DOUBLE,0,
              problem_pt->communicator_pt()->mpi_comm());
   }
#endif
  //Now construct the value
  double x2 = (psi_x1 - parameter_residual)/psi_c;

  //Now construct the vectors that multiply the jacobian terms
  Vector<DoubleVectorWithHaloEntries> D_and_X1(2);
  D_and_X1[0].build(Linear_solver_pt->distribution_pt(),0.0);
  D_and_X1[1].build(Linear_solver_pt->distribution_pt(),0.0);
  //Fill in the appropriate terms
  //Get the number of local dofs from the Linear_solver_pt distribution
  const unsigned n_dof_local_linear_solver =
   Linear_solver_pt->distribution_pt()->nrow_local();
  for(unsigned n=0;n<n_dof_local_linear_solver;n++)
   {
    const double C_ = (*C_pt)[n];
    D_and_X1[0][n] = (*D_pt)[n] - Psi*C_;
    D_and_X1[1][n] = x1[n] - x2*C_;
   }

  //Local storage for the result of the product terms
  Vector<DoubleVectorWithHaloEntries> Jprod_D_and_X1(2);

  //Redistribute the inputs to have the Dof distribution pt
  D_and_X1[0].redistribute(handler_pt->Dof_distribution_pt);
  D_and_X1[1].redistribute(handler_pt->Dof_distribution_pt);

  //Set up the halo scheme
#ifdef OOMPH_HAS_MPI
  D_and_X1[0].build_halo_scheme(problem_pt->Halo_scheme_pt);
  D_and_X1[1].build_halo_scheme(problem_pt->Halo_scheme_pt);
#endif

  //Get the products from the new problem function
  problem_pt->get_hessian_vector_products(
   handler_pt->Y,D_and_X1,Jprod_D_and_X1);

  //OK, now we can formulate the next vectors
  //(again assuming result contains residuals)
  //Need to redistribute F to the dof distribution
  //F.redistribute(handler_pt->Dof_distribution_pt);
 DoubleVector G(handler_pt->Dof_distribution_pt);

  for(unsigned n=0;n<n_dof_local;n++)
   {
    G[n] = result[n_dof_local+offset+n]
     - Jprod_D_and_X1[1][n] - x2*dRdparam[n_dof_local+offset+n];
    Jprod_D_and_X1[0][n] *= -1.0;
    Jprod_D_and_X1[0][n] -= Psi*dRdparam[n_dof_local+offset+n];
   }


  //Then redistribute back to the linear solver distribution
  G.redistribute(Linear_solver_pt->distribution_pt());
  Jprod_D_and_X1[0].redistribute(Linear_solver_pt->distribution_pt());
  Jprod_D_and_X1[1].redistribute(Linear_solver_pt->distribution_pt());

  //Linear solve to get B
  Linear_solver_pt->resolve(Jprod_D_and_X1[0],*B_pt);
  //Liner solve to get x3
  DoubleVector x3(Linear_solver_pt->distribution_pt(),0.0);
  Linear_solver_pt->resolve(G,x3);

  //Construst a couple of additional products
  double l_x3 = psi.dot(x3);
  double l_b = psi.dot(*B_pt);

  //get the last intermediate variable
  //The required parameter is that corresponding to the dof, which
  //is only stored on the root processor
  double sigma_residual = result[2*(n_dof_local+offset)-1];
#ifdef OOMPH_HAS_MPI
  //Broadcast to all others, if we have a distributed problem
  if(problem_pt->distributed())
   {
    MPI_Bcast(&sigma_residual,1,MPI_DOUBLE,0,
              problem_pt->communicator_pt()->mpi_comm());
   }
#endif


  const double delta_sigma = (l_x3 - sigma_residual)/l_b;
  const double delta_lambda = x2 - delta_sigma*Psi;

  //Need to do some more rearrangements here because result is global
  //but the other vectors are not!

  //Create temporary DoubleVectors to hold the results
  DoubleVector res1(Linear_solver_pt->distribution_pt());
  DoubleVector res2(Linear_solver_pt->distribution_pt());

  for(unsigned n=0;n<n_dof_local_linear_solver;n++)
   {
    res1[n] =
     x1[n] - delta_lambda*(*C_pt)[n] - delta_sigma*(*D_pt)[n];
    res2[n] = x3[n] - delta_sigma*(*B_pt)[n];
   }

  //Now redistribute these into the Dof distribution pointer
  res1.redistribute(handler_pt->Dof_distribution_pt);
  res2.redistribute(handler_pt->Dof_distribution_pt);

  //Now can copy over into results into the result vector
  for(unsigned n=0;n<n_dof_local;n++)
   {
    result[n] = res1[n];
    result[n_dof_local + offset + n] = res2[n];
   }

  //Add the final contributions to the residuals
  //only on the root processor if we have a distributed problem
#ifdef OOMPH_HAS_MPI
  if((!problem_pt->distributed()) || (my_rank==0))
#endif
   {
    result[n_dof_local] = delta_lambda;
    result[2*n_dof_local+1] = delta_sigma;
   }


  //The sign of the determinant is given by the sign of
  //the product psi_c and l_b
  //NOT CHECKED YET!
  problem_pt->sign_of_jacobian() =
   static_cast<int>(std::fabs(psi_c*l_b)/(psi_c*l_b));

  //Redistribute the result into its incoming distribution
  result.redistribute(&result_dist);

  //Switch things to our block solver
  handler_pt->solve_full_system();

  //If we are not storing things, clear the results
  if(!Enable_resolve)
   {
    //We no longer need to store the matrix
    Linear_solver_pt->disable_resolve();
    delete B_pt; B_pt=0;
    delete C_pt; C_pt=0;
    delete D_pt; D_pt=0;
    delete dJy_dparam_pt; dJy_dparam_pt = 0;
   }
  //Otherwise also store the pointer to the problem
  else
   {
    Problem_pt = problem_pt;
   }
 }


 //==============================================================
 //Hack the re-solve to use the block factorisation
 //==============================================================
 void BlockPitchForkLinearSolver::resolve(const DoubleVector &rhs,
                                          DoubleVector &result)
 {
     std::cout << "Block pitchfork resolve" << std::endl;
  //Check that the factors have been stored
  if(B_pt==0)
   {
    throw OomphLibError("The required vectors have not been stored",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }


  //Cache pointer to the problem
  Problem* const problem_pt = Problem_pt;

  //Locally cache pointer to the handler
  PitchForkHandler* handler_pt =
   static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());

  //Get the augmented distribution from the handler
  LinearAlgebraDistribution aug_dist =
   handler_pt->Augmented_dof_distribution_pt;
  //this->build_distribution(aug_dist);

  //Find the number of dofs of the augmented system
  //const unsigned n_aug_dof = problem_pt->ndof();

  //Storage for the result distribution
  LinearAlgebraDistribution result_dist;

  // if the result vector is not setup then rebuild with distribution
  // = natural distribution of augmented solver
  if(!result.built())
   {
    result.build(&aug_dist,0.0);
  }
  //Otherwise store the incoming distribution and redistribute
  else
   {
    result_dist.build(result.distribution_pt());
    result.redistribute(&aug_dist);
   }


  //Locally cache a pointer to the parameter
  //double* const parameter_pt = handler_pt->Parameter_pt;

  //Switch things to our block solver
  handler_pt->solve_block_system();
  //We need to find out the number of dofs
  //unsigned n_dof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  //LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
  //this->build_distribution(dist);

  // if the result vector is not setup then rebuild with distribution = global
  //if (!result.built())
  // {
  //  result.build(this->distribution_pt(),0.0);
  // }


  //Copy the rhs into local storage
  DoubleVector rhs_local = rhs;
  //and redistribute into the augmented distribution
  rhs_local.redistribute(&aug_dist);

  //Setup storage for output
  DoubleVector x1(Linear_solver_pt->distribution_pt(),0.0);
  DoubleVector x3(Linear_solver_pt->distribution_pt(),0.0);

  //Create input for RHS with the natural distribution
  DoubleVector input_x1(handler_pt->Dof_distribution_pt);
  const unsigned n_dof_local = input_x1.nrow_local();

  //Set the values of the a vector
  for(unsigned n=0;n<n_dof_local;n++) {input_x1[n] = rhs_local[n];}
  //Need to redistribute into the linear algebra distribution
  input_x1.redistribute(Linear_solver_pt->distribution_pt());

  //We want to resolve, of course
  Linear_solver_pt->enable_resolve();
  //Now solve for the first vector
  Linear_solver_pt->resolve(input_x1,x1);

  //Get the symmetry vector from the handler
  DoubleVector psi = handler_pt->Psi;
  //redistribute local copy
  psi.redistribute(Linear_solver_pt->distribution_pt());

  //We can now construct various dot products
  double psi_d = psi.dot(*D_pt);
  double psi_c = psi.dot(*C_pt);
  double psi_x1 = psi.dot(x1);

  //Calculate another intermediate constant
  const double Psi = psi_d/psi_c;

  //Construct the second intermediate value,
  //assumption is that rhs has been set to the current value of the residuals
  double parameter_residual = rhs_local[n_dof_local];
#ifdef OOMPH_HAS_MPI
  //Broadcast to all others, if we have a distributed problem
  if(problem_pt->distributed())
   {
    MPI_Bcast(&parameter_residual,1,MPI_DOUBLE,0,
              problem_pt->communicator_pt()->mpi_comm());
   }
#endif

  double x2 = (psi_x1 - parameter_residual)/psi_c;

  //Now construct the vectors that multiply the jacobian terms
  //Vector<double> X1(n_dof/*+1*/);
  Vector<DoubleVectorWithHaloEntries> X1(1);
  X1[0].build(Linear_solver_pt->distribution_pt(),0.0);

  const unsigned n_dof_local_linear_solver =
   Linear_solver_pt->distribution_pt()->nrow_local();
  for(unsigned n=0;n<n_dof_local_linear_solver;n++)
   {
    X1[0][n] = x1[n] - x2*(*C_pt)[n];
   }

  //Local storage for the product term
  Vector<DoubleVectorWithHaloEntries> Jprod_X1(1);

  //Redistribute the inputs to have the Dof distribution pt
  X1[0].redistribute(handler_pt->Dof_distribution_pt);

  //Set up the halo scheme
#ifdef OOMPH_HAS_MPI
  X1[0].build_halo_scheme(problem_pt->Halo_scheme_pt);
#endif

  //Get the product from the problem
  problem_pt->get_hessian_vector_products(handler_pt->Y,X1,Jprod_X1);

  //In standard cases the offset here will be one
  unsigned offset=1;
#ifdef OOMPH_HAS_MPI
  //If we are distributed and not on the first processor
  //then there is no offset and the eigenfunction is
  //directly after the standard dofs
  int my_rank = problem_pt->communicator_pt()->my_rank();
  if((problem_pt->distributed()) && (my_rank!=0)) {offset=0;}
#endif

  //OK, now we can formulate the next vectors
  //(again assuming result contains residuals)
  //Local storage for the product terms
  DoubleVector Mod_Jprod_X1(handler_pt->Dof_distribution_pt,0.0);

  for(unsigned n=0;n<n_dof_local;n++)
   {
    Mod_Jprod_X1[n] = rhs_local[n_dof_local+offset+n] - Jprod_X1[0][n]
     - x2*(*dJy_dparam_pt)[n];
   }

  //Redistribute back to the linear solver distribution
  Mod_Jprod_X1.redistribute(Linear_solver_pt->distribution_pt());

  //Liner solve to get x3
  Linear_solver_pt->resolve(Mod_Jprod_X1,x3);

  //Construst a couple of additional products
  double l_x3 = psi.dot(x3);
  double l_b = psi.dot(*B_pt);

  //get the last intermediate variable
  //The required parameter is that corresponding to the dof, which
  //is only stored on the root processor
  double sigma_residual = rhs_local[2*(n_dof_local+offset)-1];
#ifdef OOMPH_HAS_MPI
  //Broadcast to all others, if we have a distributed problem
  if(problem_pt->distributed())
   {
    MPI_Bcast(&sigma_residual,1,MPI_DOUBLE,0,
              problem_pt->communicator_pt()->mpi_comm());
   }
#endif

  //get the last intermediate variable
  const double delta_sigma = (l_x3 - sigma_residual)/l_b;
  const double delta_lambda = x2 - delta_sigma*Psi;

  //Create temporary DoubleVectors to hold the results
  DoubleVector res1(Linear_solver_pt->distribution_pt());
  DoubleVector res2(Linear_solver_pt->distribution_pt());

  for(unsigned n=0;n<n_dof_local_linear_solver;n++)
   {
    res1[n] = x1[n] - delta_lambda*(*C_pt)[n] - delta_sigma*(*D_pt)[n];
    res2[n] = x3[n] - delta_sigma*(*B_pt)[n];
   }

  //Now redistribute these into the Dof distribution pointer
  res1.redistribute(handler_pt->Dof_distribution_pt);
  res2.redistribute(handler_pt->Dof_distribution_pt);

  //Now can copy over into results into the result vector
  for(unsigned n=0;n<n_dof_local;n++)
   {
    result[n] = res1[n];
    result[n_dof_local + offset + n] = res2[n];
   }

  //Add the final contributions to the residuals
  //only on the root processor if we have a distributed problem
#ifdef OOMPH_HAS_MPI
  if((!problem_pt->distributed()) || (my_rank==0))
#endif
   {
    result[n_dof_local] = delta_lambda;
    result[2*n_dof_local+1] = delta_sigma;
   }

  //Redistribute the result into its incoming distribution (if it had one)
  if(result_dist.built()) {result.redistribute(&result_dist);}

  Linear_solver_pt->disable_resolve();

  //Switch things to our block solver
  handler_pt->solve_full_system();
 }


//--------------------------------------------------------------



 //======================================================================
 /// Clean up the memory that may have been allocated by the solver
 //=====================================================================
 AugmentedBlockPitchForkLinearSolver::~AugmentedBlockPitchForkLinearSolver()
 {
  if(Alpha_pt!=0) {delete Alpha_pt;}
  if(E_pt!=0) {delete E_pt;}
 }

 //===================================================================
 /// Use a block factorisation to solve the augmented system
 /// associated with a PitchFork bifurcation.
 //===================================================================
 void AugmentedBlockPitchForkLinearSolver:: solve(Problem* const &problem_pt,
                                                  DoubleVector &result)
 {
     std::cout << "Augmented pitchfork solve" << std::endl;
  // if the result is setup then it should not be distributed
#ifdef PARANOID
  if (result.built())
   {
    if (result.distributed())
     {
      throw OomphLibError("The result vector must not be distributed",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
#endif


  //Locally cache the pointer to the handler.
  PitchForkHandler* handler_pt =
   static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());

  //Switch the handler to "block solver" mode
  handler_pt->solve_augmented_block_system();

  //We need to find out the number of dofs in the problem
  unsigned n_dof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
  this->build_distribution(dist);

  // if the result vector is not setup then rebuild with distribution = global
  if (!result.built())
   {
    result.build(this->distribution_pt(),0.0);
   }

  //Setup storage for temporary vectors
  DoubleVector a(this->distribution_pt(),0.0), b(this->distribution_pt(),0.0);

  //Allocate storage for Alpha which can be used in the resolve
  if(Alpha_pt!=0) {delete Alpha_pt;}
  Alpha_pt = new DoubleVector(this->distribution_pt(),0.0);

  //We are going to do resolves using the underlying linear solver
  Linear_solver_pt->enable_resolve();
  //Solve the first system Aa = R
  Linear_solver_pt->solve(problem_pt,a);

  //Get the symmetry vector from the handler
  DoubleVector psi(this->distribution_pt(),0.0);
  for(unsigned n=0;n<(n_dof-1);++n)  {psi[n] = handler_pt->Psi[n];}
  //Set the final entry to zero
  psi[n_dof-1] = 0.0;

  //Now resolve to find alpha
  Linear_solver_pt->resolve(psi,*Alpha_pt);

  //We can now construct our multipliers
  //Prepare to scale
  double dof_length=0.0, a_length=0.0, alpha_length=0.0;
  for(unsigned n=0;n<n_dof;n++)
   {
    if(std::fabs(problem_pt->dof(n)) > dof_length)
     {dof_length = std::fabs(problem_pt->dof(n));}
    if(std::fabs(a[n]) > a_length) {a_length = std::fabs(a[n]);}
    if(std::fabs((*Alpha_pt)[n]) > alpha_length)
     {alpha_length = std::fabs((*Alpha_pt)[n]);}
   }

  double a_mult = dof_length/a_length;
  double alpha_mult = dof_length/alpha_length;
  const double FD_step = 1.0e-8;
  a_mult += FD_step; alpha_mult += FD_step;
  a_mult *= FD_step; alpha_mult *= FD_step;

  //Local storage for the product terms
  DoubleVector Jprod_a(this->distribution_pt(),0.0),
   Jprod_alpha(this->distribution_pt(),0.0);

  //Calculate the product of the jacobian matrices, etc
  unsigned long n_element = problem_pt->mesh_pt()->nelement();
  for(unsigned long e = 0;e<n_element;e++)
   {
    GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the ndofs in each element
    unsigned n_var = handler_pt->ndof(elem_pt);
    //Get the jacobian matrices
    DenseMatrix<double> jac(n_var), jac_a(n_var), jac_alpha(n_var);
     // the elemental residual
     Vector<double> res_elemental(n_var);
    //Get unperturbed jacobian
    handler_pt->get_jacobian(elem_pt,res_elemental,jac);

    //Backup the dofs
    Vector<double> dof_bac(n_var);
    //Perturb the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
      dof_bac[n] = problem_pt->dof(eqn_number);
      //Pertub by vector a
      problem_pt->dof(eqn_number) += a_mult*a[eqn_number];
     }

    problem_pt->actions_after_change_in_bifurcation_parameter();

    //Now get the new jacobian
    handler_pt->get_jacobian(elem_pt,res_elemental,jac_a);

    //Perturb the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
      problem_pt->dof(eqn_number) = dof_bac[n];
      //Pertub by vector a
      problem_pt->dof(eqn_number) += alpha_mult*(*Alpha_pt)[eqn_number];
     }

    problem_pt->actions_after_change_in_bifurcation_parameter();

    //Now get the new jacobian
    handler_pt->get_jacobian(elem_pt,res_elemental,jac_alpha);

    //Reset the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
      problem_pt->dof(eqn_number) = dof_bac[n];
     }

    problem_pt->actions_after_change_in_bifurcation_parameter();

    //OK, now work out the products
    for(unsigned n=0;n<(n_var-1);n++)
     {
      unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
      double prod_a=0.0, prod_alpha=0.0;
      for(unsigned m=0;m<(n_var-1);m++)
       {
        unsigned unknown = handler_pt->eqn_number(elem_pt,m);
        prod_a += (jac_a(n,m) - jac(n,m))*
         handler_pt->Y[unknown];
        prod_alpha += (jac_alpha(n,m) - jac(n,m))*
         handler_pt->Y[unknown];
       }
      Jprod_a[eqn_number] += prod_a/a_mult;
      Jprod_alpha[eqn_number] += prod_alpha/alpha_mult;
     }
   }

   Jprod_alpha[n_dof-1] = 0.0;
   Jprod_a[n_dof-1] = 0.0;

   //OK, now we can formulate the next vectors
   //The assumption here is that the result has been set to the
   //residuals.
   for(unsigned n=0;n<n_dof-1;n++)
    {
     b[n] = result[n_dof+n] - Jprod_a[n];
    }
   b[n_dof-1] = result[2*n_dof-1];


   //Allocate storage for E which can be used in the resolve
   if(E_pt!=0) {delete E_pt;}
   E_pt = new DoubleVector(this->distribution_pt(),0.0);

   DoubleVector f(this->distribution_pt(),0.0);

   Linear_solver_pt->resolve(b,f);
   Linear_solver_pt->resolve(Jprod_alpha,*E_pt);

   //Calculate the final entry in the vector e
   const double e_final = (*E_pt)[n_dof-1];
   //Calculate the final entry in the vector d
   const double d_final = -f[n_dof-1]/e_final;
   //Assemble the final corrections
   for(unsigned n=0;n<n_dof-1;n++)
    {
     result[n] = a[n] - (*Alpha_pt)[n]*d_final;
     result[n_dof+n] = f[n] + (*E_pt)[n]*d_final;
    }

   result[n_dof-1] = a[n_dof-1] - (*Alpha_pt)[n_dof-1]*d_final;
   result[2*n_dof-1] = d_final;

   //The sign of the jacobian is the sign of the final entry in e
   problem_pt->sign_of_jacobian() =
    static_cast<int>(std::fabs(e_final)/e_final);


   //Switch things to our block solver
   handler_pt->solve_full_system();

   //If we are not storing things, clear the results
   if(!Enable_resolve)
    {
     //We no longer need to store the matrix
     Linear_solver_pt->disable_resolve();
     delete Alpha_pt; Alpha_pt=0;
     delete E_pt; E_pt=0;
    }
   //Otherwise also store the pointer to the problem
   else
    {
     Problem_pt = problem_pt;
    }
 }


 //==============================================================
 //Hack the re-solve to use the block factorisation
 //==============================================================
 void AugmentedBlockPitchForkLinearSolver::resolve(const DoubleVector & rhs,
                                                   DoubleVector &result)
 {
     std::cout << "Augmented pitchfork resolve" << std::endl;
  //Check that the factors have been stored
  if(Alpha_pt==0)
   {
    throw OomphLibError("The required vectors have not been stored",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

  //Get the pointer to the problem
  Problem* const problem_pt = Problem_pt;

  PitchForkHandler* handler_pt =
   static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());

  //Switch things to our block solver
  handler_pt->solve_augmented_block_system();
  //We need to find out the number of dofs
  unsigned n_dof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
  this->build_distribution(dist);

  // if the result vector is not setup then rebuild with distribution = global
  if (!result.built())
   {
    result.build(this->distribution_pt(),0.0);
   }


  //Setup storage
  DoubleVector a(this->distribution_pt(),0.0), b(this->distribution_pt(),0.0);

  //Set the values of the a vector
  for(unsigned n=0;n<n_dof;n++) {a[n] = rhs[n];}

  Linear_solver_pt->enable_resolve();

  // Copy rhs vector into local storage so it doesn't get overwritten
  // if the linear solver decides to initialise the solution vector, say,
  // which it's quite entitled to do!
  DoubleVector input_a(a);

  Linear_solver_pt->resolve(input_a,a);

  //We can now construct our multipliers
  //Prepare to scale
  double dof_length=0.0, a_length=0.0;
  for(unsigned n=0;n<n_dof;n++)
   {
    if(std::fabs(problem_pt->dof(n)) > dof_length)
     {dof_length = std::fabs(problem_pt->dof(n));}

    if(std::fabs(a[n]) > a_length) {a_length = std::fabs(a[n]);}
   }
  double a_mult = dof_length/a_length;
  const double FD_step = 1.0e-8;
  a_mult += FD_step;
  a_mult *= FD_step;

  DoubleVector Jprod_a(this->distribution_pt(),0.0);

  unsigned long n_element = problem_pt->mesh_pt()->nelement();
   for(unsigned long e = 0;e<n_element;e++)
    {
     GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
     //Loop over the ndofs in each element
     unsigned n_var = handler_pt->ndof(elem_pt);
     //Get some jacobian matrices
     DenseMatrix<double> jac(n_var), jac_a(n_var);
     // the elemental residual
     Vector<double> res_elemental(n_var);
     //Get unperturbed jacobian
     handler_pt->get_jacobian(elem_pt,res_elemental,jac);

     //Backup the dofs
     DoubleVector dof_bac(this->distribution_pt(),0.0);
     //Perturb the dofs
     for(unsigned n=0;n<n_var;n++)
      {
       unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
       dof_bac[n] = problem_pt->dof(eqn_number);
       //Pertub by vector a
       problem_pt->dof(eqn_number) += a_mult*a[eqn_number];
      }

     problem_pt->actions_after_change_in_bifurcation_parameter();

     //Now get the new jacobian
     handler_pt->get_jacobian(elem_pt,res_elemental,jac_a);

     //Reset the dofs
     for(unsigned n=0;n<n_var;n++)
      {
       unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
       problem_pt->dof(eqn_number) = dof_bac[n];
      }

     problem_pt->actions_after_change_in_bifurcation_parameter();

     //OK, now work out the products
     for(unsigned n=0;n<(n_var-1);n++)
      {
       unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
       double prod_a=0.0;
       for(unsigned m=0;m<(n_var-1);m++)
        {
         unsigned unknown = handler_pt->eqn_number(elem_pt,m);
         prod_a += (jac_a(n,m) - jac(n,m))*
          handler_pt->Y[unknown];
        }
       Jprod_a[eqn_number] += prod_a/a_mult;
      }
    }

   Jprod_a[n_dof-1] = 0.0;

   //OK, now we can formulate the next vectors
   for(unsigned n=0;n<n_dof-1;n++)
    {
     b[n] = rhs[n_dof+n] - Jprod_a[n];
    }
   b[n_dof-1] = rhs[2*n_dof-1];

   DoubleVector f(this->distribution_pt(),0.0);

   Linear_solver_pt->resolve(b,f);

   //Calculate the final entry in the vector d
   const double d_final = -f[n_dof-1]/(*E_pt)[n_dof-1];
   //Assemble the final corrections
   for(unsigned n=0;n<n_dof-1;n++)
    {
     result[n] = a[n] - (*Alpha_pt)[n]*d_final;
     result[n_dof+n] = f[n] + (*E_pt)[n]*d_final;
    }

   result[n_dof-1] = a[n_dof-1] - (*Alpha_pt)[n_dof-1]*d_final;
   result[2*n_dof-1] = d_final;

   Linear_solver_pt->disable_resolve();

   //Switch things to our block solver
   handler_pt->solve_full_system();
  }

  ///////////////////////////////////////////////////////////////////////
 // Non-inline functions for the PitchForkHandler class
 //////////////////////////////////////////////////////////////////////

 //==================================================================
 ///Constructor: Initialise the PitchForkHandler by setting intial
 ///guesses for Sigma, Y, specifying C and Psi and calculating count.
 //==================================================================
 PitchForkHandler::PitchForkHandler(Problem* const &problem_pt,
                                    AssemblyHandler* const
                                    &assembly_handler_pt,
                                    double* const &parameter_pt,
                                    const DoubleVector &symmetry_vector) :
  Solve_which_system(Full_augmented),
  Sigma(0.0),
  Parameter_pt(parameter_pt)
 {
  //Set the problem pointer
  Problem_pt = problem_pt;
  //Set the assembly handler
  Assembly_handler_pt = assembly_handler_pt;
  //Set the number of degrees of freedom
  Ndof = Problem_pt->ndof();
  //Backup the distribution
  Dof_distribution_pt  = Problem_pt->Dof_distribution_pt;
#ifdef OOMPH_HAS_MPI
  //Set the distribution flag
  Distributed = Problem_pt->distributed();
#endif

  //Find the elements in the problem
  unsigned n_element = Problem_pt->mesh_pt()->nelement();

#ifdef OOMPH_HAS_MPI

  //Work out all the global equations to which this processor
  //contributes and store the halo data in the problem
  Problem_pt->setup_dof_halo_scheme();

#endif


  //Now use the dof distribution for all double vectors
  Psi.build(Dof_distribution_pt);
  Y.build(Dof_distribution_pt);
  C.build(Dof_distribution_pt);
  Count.build(Dof_distribution_pt);
#ifdef OOMPH_HAS_MPI
  Count.build_halo_scheme(Problem_pt->Halo_scheme_pt);
#endif

  //Loop over the elements and count the entries
  //and number of (non-halo) elements
  unsigned n_non_halo_element_local=0;
  for(unsigned e=0;e<n_element;e++)
   {
    GeneralisedElement* elem_pt = Problem_pt->mesh_pt()->element_pt(e);
#ifdef OOMPH_HAS_MPI
    //Ignore halo elements
    if(!elem_pt->is_halo())
     {
#endif
      //Increment the number of non halo elements
      ++n_non_halo_element_local;
      //Now count the number of times the element contributes to a value
      unsigned n_var = assembly_handler_pt->ndof(elem_pt);
      for(unsigned n=0;n<n_var;n++)
       {
        ++Count.global_value(assembly_handler_pt->eqn_number(elem_pt,n));
       }
#ifdef OOMPH_HAS_MPI
     }
#endif
   }

  //Add together all the counts
#ifdef OOMPH_HAS_MPI
  Count.sum_all_halo_and_haloed_values();

  //If distributed, find the total number of elements in the problem
  if(Distributed)
   {
    //Need to gather the total number of non halo elements
    MPI_Allreduce(&n_non_halo_element_local,&Nelement,1,MPI_UNSIGNED,MPI_SUM,
                  Problem_pt->communicator_pt()->mpi_comm());
   }
  //Otherwise the total number is the same on each processor
  else
#endif
   {
    Nelement = n_non_halo_element_local;
   }

#ifdef OOMPH_HAS_MPI
  //Only add the parameter to the first processor, if distributed
  int my_rank = Problem_pt->communicator_pt()->my_rank();
  if((!Distributed) || (my_rank==0))
#endif
   {
    //Add the parameter to the problem
    Problem_pt->Dof_pt.push_back(parameter_pt);
   }

  //Find length of the symmetry vector
  double length = symmetry_vector.norm();

  //Add the unknowns for the null vector to the problem
  //Normalise the symmetry vector and initialise the null vector to the
  //symmetry vector and set the constant vector, c, to the
  //normalised symmetry vector.

  //Need to redistribute the symmetry vector into the (natural)
  //distribution of the Dofs
#ifdef OOMPH_HAS_MPI
  //Check that the symmetry vector has the correct distribution
  if(!(*symmetry_vector.distribution_pt() == *Dof_distribution_pt))
   {
    throw OomphLibError
     ("The symmetry vector must have the same distribution as the dofs\n",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
   }
#endif

  //Only loop over the local unknowns
  const unsigned n_dof_local = Dof_distribution_pt->nrow_local();
  for(unsigned n=0;n<n_dof_local;n++)
   {
    Problem_pt->Dof_pt.push_back(&Y[n]);
    //Psi[n] = symmetry_vector[n];
    Psi[n] = Y[n] = C[n] = symmetry_vector[n]/length;
   }


#ifdef OOMPH_HAS_MPI
  //Set up the required halo schemes (which also synchronises)
  Psi.build_halo_scheme(Problem_pt->Halo_scheme_pt);
  Y.build_halo_scheme(Problem_pt->Halo_scheme_pt);
  C.build_halo_scheme(Problem_pt->Halo_scheme_pt);


  if((!Distributed) || (my_rank==0))
#endif
   //Add the slack parameter to the problem on the first processor
   //if distributed
   {
    Problem_pt->Dof_pt.push_back(&Sigma);
   }

#ifdef OOMPH_HAS_MPI
  if(Distributed)
   {
    unsigned augmented_first_row=0;
    unsigned augmented_n_row_local=0;

    //Set up the translation scheme on every processor
    Global_eqn_number.resize(2*Ndof+2);
    int n_proc = Problem_pt->communicator_pt()->nproc();
    unsigned global_eqn_count=0;
    for(int d=0;d<n_proc;d++)
     {
      //Find out the first row of the current processor
      if(my_rank==d) {augmented_first_row = global_eqn_count;}

      const unsigned n_row_local = Dof_distribution_pt->nrow_local(d);
      const unsigned first_row = Dof_distribution_pt->first_row(d);
      //Add the basic equations
      for(unsigned n=0;n<n_row_local;n++)
       {
        Global_eqn_number[first_row+n] = global_eqn_count;
        ++global_eqn_count;
       }
      //If on the first processor add the pointer to the parameter
      if(d==0)
       {
        Global_eqn_number[Ndof] = global_eqn_count;
        ++global_eqn_count;
       }
      //Add the eigenfunction
      for(unsigned n=0;n<n_row_local;n++)
       {
        Global_eqn_number[Ndof + 1 + first_row+n] = global_eqn_count;
        ++global_eqn_count;
       }
      //Finally add the slack parameter
      if(d==0)
       {
        Global_eqn_number[2*Ndof+1] = global_eqn_count;
        ++global_eqn_count;
       }
      //Find out the number of rows of the current processor
      if(my_rank==d)
       {augmented_n_row_local = global_eqn_count - augmented_first_row;}
     }

    //Make a new linear algebra distribution
    Augmented_dof_distribution_pt =
     new LinearAlgebraDistribution(Problem_pt->communicator_pt(),
                                   augmented_first_row,
                                   augmented_n_row_local);
   }
  else
#endif
   {
    Augmented_dof_distribution_pt = new
     LinearAlgebraDistribution(Problem_pt->communicator_pt(),
                               2*Ndof+2,false);
   }

    // resize
    //Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
    //                                       2*Ndof+2,false);

  Problem_pt->Dof_distribution_pt = Augmented_dof_distribution_pt;

  //Remove all previous sparse storage used during Jacobian assembly
  Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 }

 //=====================================================================
 /// Destructor return the problem to its original (non-augmented) state
 //=====================================================================
 PitchForkHandler::~PitchForkHandler()
 {
  //If we are using the block solver reset the problem's linear solver
  //to the original one
  BlockPitchForkLinearSolver* block_pitchfork_solver_pt
   = dynamic_cast<BlockPitchForkLinearSolver*>(
    Problem_pt->linear_solver_pt());

  if(block_pitchfork_solver_pt)
   {
    //Reset the problem's linear solver
    Problem_pt->linear_solver_pt() =
     block_pitchfork_solver_pt->linear_solver_pt();
    //Delete the block solver
    delete block_pitchfork_solver_pt;
   }

  //If we are using the augmented
  //block solver reset the problem's linear solver
  //to the original one
  AugmentedBlockPitchForkLinearSolver* augmented_block_pitchfork_solver_pt
   = dynamic_cast<AugmentedBlockPitchForkLinearSolver*>(
    Problem_pt->linear_solver_pt());

  if(augmented_block_pitchfork_solver_pt)
   {
    //Reset the problem's linear solver
    Problem_pt->linear_solver_pt() =
     augmented_block_pitchfork_solver_pt->linear_solver_pt();
    //Delete the block solver
    delete augmented_block_pitchfork_solver_pt;
   }

  //Now return the problem to its original size
  Problem_pt->Dof_pt.resize(this->Dof_distribution_pt->nrow_local());

  //Problem_pt->Dof_pt.resize(Ndof);
  //Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
  //                                       Ndof,false);

  Problem_pt->Dof_distribution_pt = this->Dof_distribution_pt;
  //Remove all previous sparse storage used during Jacobian assembly
  Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 }

 //================================================================
 ///Get the number of elemental degrees of freedom
 //================================================================
 unsigned PitchForkHandler::ndof(GeneralisedElement* const &elem_pt)
 {
  unsigned raw_ndof = elem_pt->ndof();

  //Return different values depending on the type of block decomposition
  switch(Solve_which_system)
   {
   case Full_augmented:
    return (2*raw_ndof + 2);
    break;

   case Block_augmented_J:
    return (raw_ndof+1);
    break;

   case Block_J:
    return raw_ndof;
    break;

   default:
    std::ostringstream error_stream;
    error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
                 << " not " << Solve_which_system << "\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 }

 ///Get the global equation number of the local unknown
 unsigned long PitchForkHandler::eqn_number(GeneralisedElement* const &elem_pt,
                                            const unsigned &ieqn_local)
 {
  //Get the raw value
  unsigned raw_ndof = elem_pt->ndof();
  unsigned long global_eqn=0;

  //Which system are we solving
  switch(Solve_which_system)
   {
    //In the block case, it's just the standard numbering
   case Block_J:
    global_eqn = elem_pt->eqn_number(ieqn_local);
    break;

    //Block augmented not done properly yet
   case Block_augmented_J:
    //Not done the distributed case
#ifdef OOMPH_HAS_MPI
    if(Distributed)
     {
      throw OomphLibError(
       "Block Augmented solver not implemented for distributed case\n",
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
#endif

    //Full case
   case Full_augmented:
    //The usual equations
    if(ieqn_local < raw_ndof)
     {
      global_eqn = this->global_eqn_number(elem_pt->eqn_number(ieqn_local));
     }
    //The bifurcation parameter equation
    else if(ieqn_local == raw_ndof)
     {
      global_eqn = this->global_eqn_number(Ndof);
     }
    //If we are assembling the full system we also have
    //The components of the null vector
    else if(ieqn_local < (2*raw_ndof + 1))
     {
      global_eqn = this->global_eqn_number(
       Ndof + 1 + elem_pt->eqn_number(ieqn_local - 1 - raw_ndof));
     }
    //The slack parameter
    else
     {
      global_eqn = this->global_eqn_number(2*Ndof+1);
     }
    break;
   } //End of switch
  return global_eqn;
 }

 //==============================================================
 ///Get the residuals
 //==============================================================
 void PitchForkHandler::get_residuals(GeneralisedElement* const &elem_pt,
                                      Vector<double> &residuals)
 {

  //Need to get raw residuals and jacobian
  unsigned raw_ndof = elem_pt->ndof();

  //Find out which system we are solving
  switch(Solve_which_system)
   {
    //If we are solving the original system
   case Block_J:
   {
    //get the basic residuals
    elem_pt->get_residuals(residuals);
    //Now multiply to fill in the residuals for the final term
    for(unsigned i=0;i<raw_ndof;i++)
     {
      unsigned local_eqn = elem_pt->eqn_number(i);
      //Add the slack parameter to the final residuals
      residuals[i] += Sigma*Psi.global_value(local_eqn)/
       Count.global_value(local_eqn);
     }
   }
   break;

   //If we are solving the augmented-by-one system
   case Block_augmented_J:
   {
    //Get the basic residuals
    elem_pt->get_residuals(residuals);

    //Zero the final residual
    residuals[raw_ndof] = 0.0;
    //Now multiply to fill in the residuals for the final term
    for(unsigned i=0;i<raw_ndof;i++)
     {
      unsigned local_eqn = elem_pt->eqn_number(i);
      //Add the slack parameter to the final residuals
      residuals[i] += Sigma*Psi.global_value(local_eqn)/
       Count.global_value(local_eqn);
      //Final term that specifies the symmetry
      residuals[raw_ndof] +=
       ((*Problem_pt->global_dof_pt(local_eqn))*Psi.global_value(local_eqn))/
       Count.global_value(local_eqn);
     }
   }
   break;

   //Otherwise we are solving the fully augemented system
   case Full_augmented:
   {
    DenseMatrix<double> jacobian(raw_ndof);

    //Get the basic residuals and jacobian
    elem_pt->get_jacobian(residuals,jacobian);

    //Initialise the final residuals
    residuals[raw_ndof] = 0.0;
    residuals[2*raw_ndof+1] = -1.0/Nelement;

    //Now multiply to fill in the residuals associated
    //with the null vector condition
    for(unsigned i=0;i<raw_ndof;i++)
     {
      unsigned local_eqn = elem_pt->eqn_number(i);
      residuals[raw_ndof+1+i] = 0.0;
      for(unsigned j=0;j<raw_ndof;j++)
       {
        unsigned local_unknown = elem_pt->eqn_number(j);
        residuals[raw_ndof+1+i] += jacobian(i,j)*Y.global_value(local_unknown);
       }

      //Add the slack parameter to the governing equations
      residuals[i] += Sigma*Psi.global_value(local_eqn)/
       Count.global_value(local_eqn);

      //Specify the symmetry
      residuals[raw_ndof] +=
       ((*Problem_pt->global_dof_pt(local_eqn))*Psi.global_value(local_eqn))/
       Count.global_value(local_eqn);
      //Specify the normalisation
      residuals[2*raw_ndof+1] += (Y.global_value(local_eqn)*
                                  C.global_value(local_eqn))/
       Count.global_value(local_eqn);
     }
   }
   break;

   default:
    std::ostringstream error_stream;
    error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
                 << " not " << Solve_which_system << "\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

 }

  //======================================================================
  /// \short Calculate the elemental Jacobian matrix "d equation
  /// / d variable".
  //======================================================================
  void PitchForkHandler::get_jacobian(GeneralisedElement* const &elem_pt,
                    Vector<double> &residuals,
                    DenseMatrix<double> &jacobian)
   {

    unsigned augmented_ndof = ndof(elem_pt);
    unsigned raw_ndof = elem_pt->ndof();

    //Which system are we solving
    switch(Solve_which_system)
     {
      //If we are solving the original system
     case Block_J:
     {
      //get the raw residuals and jacobian
      elem_pt->get_jacobian(residuals,jacobian);

      //Now multiply to fill in the residuals for the final term
      for(unsigned i=0;i<raw_ndof;i++)
       {
        unsigned local_eqn = elem_pt->eqn_number(i);
        //Add the slack parameter to the final residuals
        residuals[i] += Sigma*Psi.global_value(local_eqn)/
         Count.global_value(local_eqn);
       }
     }
     break;

     //If we are solving the augmented-by-one system
     case Block_augmented_J:
     {
      //Get the full residuals, we need them
      get_residuals(elem_pt,residuals);

      //Need to get the raw jacobian (and raw residuals)
      Vector<double> newres(augmented_ndof);
      elem_pt->get_jacobian(newres,jacobian);

      //Now do finite differencing stuff
      const double FD_step = 1.0e-8;
      //Fill in the first lot of finite differences
      {
       //increase the global parameter
       double *unknown_pt = Parameter_pt;
       double init = *unknown_pt;
       *unknown_pt += FD_step;

       //Not do any possible updates
       Problem_pt
        ->actions_after_change_in_bifurcation_parameter();

       //Get the new (modified) residuals
       get_residuals(elem_pt,newres);

       //The final column  is given by the difference
       //between the residuals
       for(unsigned n=0;n<raw_ndof;n++)
        {
         jacobian(n,augmented_ndof-1) = (newres[n] - residuals[n])/FD_step;
        }
       //Reset the global parameter
       *unknown_pt = init;

       //Now do any possible updates
       Problem_pt
        ->actions_after_change_in_bifurcation_parameter();
      }

      //Fill in the bottom row
      for(unsigned n=0;n<raw_ndof;n++)
       {
        unsigned local_eqn = elem_pt->eqn_number(n);
        jacobian(augmented_ndof-1,n)
         = Psi.global_value(local_eqn)/Count.global_value(local_eqn);
       }
     }
     break;

     //Otherwise solving the full system
     case Full_augmented:
     {
         /// ALICE:
         Problem_pt->actions_before_newton_convergence_check();
      //Get the basic jacobian and residuals
      elem_pt->get_jacobian(residuals,jacobian);
      //get the proper residuals
      get_residuals(elem_pt,residuals);

      //Now fill in the next diagonal jacobian entry
      for(unsigned n=0;n<raw_ndof;n++)
       {
        for(unsigned m=0;m<raw_ndof;m++)
         {
          jacobian(raw_ndof+1+n,raw_ndof+1+m) = jacobian(n,m);
         }
        unsigned local_eqn = elem_pt->eqn_number(n);
        //Add in the sigma contribution
        jacobian(n,2*raw_ndof+1) = Psi.global_value(local_eqn)/
         Count.global_value(local_eqn);
        //Symmetry constraint
        jacobian(raw_ndof,n) = Psi.global_value(local_eqn)/
         Count.global_value(local_eqn);
        //Non-zero constraint
        jacobian(2*raw_ndof+1,raw_ndof+1+n) = C.global_value(local_eqn)/
         Count.global_value(local_eqn);
       }

      //Finite difference the remaining blocks
      const double FD_step = 1.0e-8;

      Vector<double> newres_p(augmented_ndof);

      //Loop over the ndofs only
      for(unsigned n=0;n<raw_ndof;++n)
       {
        //Get the original (unaugmented) global equation number
        unsigned long global_eqn = Assembly_handler_pt->eqn_number(elem_pt,n);
        double* unknown_pt = Problem_pt->global_dof_pt(global_eqn);
        double init = *unknown_pt;
        *unknown_pt += FD_step;
Problem_pt->actions_before_newton_convergence_check();
        //Get the new residuals
        get_residuals(elem_pt,newres_p);
        //Fill in the entries in the block d(Jy)/dx
        for(unsigned m=0;m<raw_ndof;m++)
         {
          jacobian(raw_ndof+1+m,n) =
           (newres_p[raw_ndof+1+m] - residuals[raw_ndof+1+m])/(FD_step);
         }

        //Reset the unknown
        *unknown_pt = init;
        Problem_pt->actions_before_newton_convergence_check();
       }

      {
       //Now increase the global parameter
       double* unknown_pt = Parameter_pt;
       double init = *unknown_pt;
       *unknown_pt += FD_step;

       Problem_pt->actions_after_change_in_bifurcation_parameter();
       //Get the new residuals
       get_residuals(elem_pt,newres_p);

       //Add in the first block d/dx
       for(unsigned m=0;m<raw_ndof;m++)
        {
         jacobian(m,raw_ndof) =
            (newres_p[m] - residuals[m])/FD_step;
        }
       //Add in the second block d/dy
       for(unsigned m=raw_ndof+1;m<augmented_ndof-1;m++)
        {
         jacobian(m,raw_ndof) =
          (newres_p[m] - residuals[m])/FD_step;
        }

       //Reset the unknown
       *unknown_pt = init;
       Problem_pt->actions_after_change_in_bifurcation_parameter();
      }
     }
     break;

     default:
      std::ostringstream error_stream;
      error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
                   << " not " << Solve_which_system << "\n";
      throw OomphLibError(error_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
        /// ALICE
              Problem_pt
        ->actions_after_change_in_bifurcation_parameter();
   }



 //==============================================================
 ///Get the derivatives of the residuals with respect to a parameter
 //==============================================================
 void PitchForkHandler::get_dresiduals_dparameter(
  GeneralisedElement* const &elem_pt, double* const &parameter_pt,
  Vector<double> &dres_dparam)
 {
  //Need to get raw residuals and jacobian
  unsigned raw_ndof = elem_pt->ndof();
Problem_pt->actions_before_newton_convergence_check();
  //Find out which system we are solving
  switch(Solve_which_system)
   {
    //If we are solving the original system
   case Block_J:
   {
    //get the basic residual derivatives
    elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);
    //Slack parameter term does not depened explicitly on the parameter
    //so is not added
   }
   break;

   //If we are solving the augmented-by-one system
   case Block_augmented_J:
   {
    //Get the basic residual derivatives
    elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);

    //Zero the final residuals derivative
    dres_dparam[raw_ndof] = 0.0;

    //Other terms must not depend on the parameter
   }
   break;

   //Otherwise we are solving the fully augemented system
   case Full_augmented:
   {
    DenseMatrix<double> djac_dparam(raw_ndof);

    //Get the basic residuals and jacobian derivatives
    elem_pt->get_djacobian_dparameter(parameter_pt,dres_dparam,djac_dparam);

    //The "final" residual derivatives do not depend on the parameter
    dres_dparam[raw_ndof] = 0.0;
    dres_dparam[2*raw_ndof+1] = 0.0;

    //Now multiply to fill in the residuals associated
    //with the null vector condition
    for(unsigned i=0;i<raw_ndof;i++)
     {
      dres_dparam[raw_ndof+1+i] = 0.0;
      for(unsigned j=0;j<raw_ndof;j++)
       {
        unsigned local_unknown = elem_pt->eqn_number(j);
        dres_dparam[raw_ndof+1+i] +=
         djac_dparam(i,j)*Y.global_value(local_unknown);
       }
     }
   }
   break;

   default:
    std::ostringstream error_stream;
    error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
                 << " not " << Solve_which_system << "\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

 }



 //========================================================================
 /// Overload the derivative of the residuals and jacobian
 /// with respect to a parameter so that it breaks because it should not
 /// be required
 //========================================================================
 void PitchForkHandler::get_djacobian_dparameter(
  GeneralisedElement* const &elem_pt,
  double* const &parameter_pt,
  Vector<double> &dres_dparam,
  DenseMatrix<double> &djac_dparam)
 {
  std::ostringstream error_stream;
  error_stream <<
   "This function has not been implemented because it is not required\n";
  error_stream << "in standard problems.\n";
  error_stream <<
   "If you find that you need it, you will have to implement it!\n\n";

  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
 }


 //=====================================================================
 /// Overload the hessian vector product function so that
 /// it calls the underlying element's hessian function.
 //========================================================================
 void PitchForkHandler::get_hessian_vector_products(
  GeneralisedElement* const &elem_pt,
  Vector<double> const &Y,
  DenseMatrix<double> const &C,
  DenseMatrix<double> &product)
 {
  elem_pt->get_hessian_vector_products(Y,C,product);
 }



 //==========================================================================
 /// Return the eigenfunction(s) associated with the bifurcation that
 /// has been detected in bifurcation tracking problems
 //==========================================================================
 void PitchForkHandler::get_eigenfunction(
  Vector<DoubleVector> &eigenfunction)
 {
  //There is only one (real) null vector
  eigenfunction.resize(1);
  //Rebuild the vector
  eigenfunction[0].build(this->Dof_distribution_pt,0.0);
  //Set the value to be the null vector
  const unsigned n_row_local = eigenfunction[0].nrow_local();
  for(unsigned n=0;n<n_row_local;n++)
   {
    eigenfunction[0][n] = Y[n];
   }
 }

#ifdef OOMPH_HAS_MPI
//=====================================================================
//Synchronise the required data
//====================================================================
 void  PitchForkHandler::synchronise()
 {
  //Only need to bother if the problem is distributed
  if(Distributed)
   {
    //Need only to synchronise the eigenfunction
    Y.synchronise();
    //Also need to synchronise the parameter and the slack parameter
    double broadcast_data[2];
    broadcast_data[0] = *Parameter_pt;
    broadcast_data[1] = Sigma;
    //Broadcast from the root to all processors
    MPI_Bcast(broadcast_data,2,MPI_DOUBLE,0,
              Problem_pt->communicator_pt()->mpi_comm());

    //Now copy received the values back
    *Parameter_pt = broadcast_data[0];
    Sigma = broadcast_data[1];
   }
 }
#endif

 //====================================================================
 /// Set to solve the augmented-by-one block system.
 //===================================================================
 void PitchForkHandler::solve_augmented_block_system()
 {
  //Only bother to do anything if we haven't already set the flag
  if(Solve_which_system!=Block_augmented_J)
   {
    //If we were solving the system with the original jacobian add the
    //parameter back
    if(Solve_which_system==Block_J)
     {Problem_pt->Dof_pt.push_back(Parameter_pt);}

    //Restrict the problem to the standard variables and
    //the bifurcation parameter only
    Problem_pt->Dof_distribution_pt = new LinearAlgebraDistribution
     (Problem_pt->communicator_pt(),Ndof+1,false);
    Problem_pt->Dof_pt.resize(Ndof+1);

//Problem_pt->Dof_pt.resize(Ndof+1);
//    Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
//                                         Ndof+1,false);

    //Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);

    //Set the solve flag
    Solve_which_system = Block_augmented_J;
   }
 }



  //==============================================================
  /// Set to solve the block system. The boolean flag is used
  /// to specify whether the block decomposition should use exactly
  /// the same jacobian as the original system.
  //==============================================================
  void PitchForkHandler::solve_block_system()
   {
    //Only bother to do anything if we haven't already set the flag
    if(Solve_which_system!=Block_J)
     {
      //Restrict the problem to the standard variables
      Problem_pt->Dof_pt.resize(this->Dof_distribution_pt->nrow_local());
      Problem_pt->Dof_distribution_pt = this->Dof_distribution_pt;

      //Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
      //Ndof,false);

      //Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);

      //Set the solve flag
      Solve_which_system = Block_J;
     }
   }

  //==============================================================
  /// \short Solve non-block system
  //==============================================================
 void PitchForkHandler::solve_full_system()
   {
    //Only do something if we are not solving the full system
    if(Solve_which_system!=Full_augmented)
     {
#ifdef OOMPH_HAS_MPI
      int my_rank = Problem_pt->communicator_pt()->my_rank();
#endif
      //If we were solving the problem without any augementation
      //add the parameter again
      if(Solve_which_system==Block_J)
       {
#ifdef OOMPH_HAS_MPI

        if((!Distributed) || (my_rank==0))
#endif
         {
          Problem_pt->Dof_pt.push_back(Parameter_pt);
         }
       }

      //Always add the additional unknowns back into the problem
      const unsigned n_dof_local = Dof_distribution_pt->nrow_local();
      for(unsigned n=0;n<n_dof_local;n++)
       {
        Problem_pt->Dof_pt.push_back(&Y[n]);
       }



#ifdef OOMPH_HAS_MPI
      if((!Distributed) || (my_rank==0))
#endif
       //Add the slack parameter to the problem on the first processor
       //if distributed
       {
        Problem_pt->Dof_pt.push_back(&Sigma);
       }


      //Delete the distribtion created in the augmented block solve
      if(Problem_pt->Dof_distribution_pt!=this->Dof_distribution_pt)
       {delete Problem_pt->Dof_distribution_pt;}

      // update the Dof distribution
      Problem_pt->Dof_distribution_pt =
       this->Augmented_dof_distribution_pt;
      //Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
      //                                       Ndof*2+2,false);
      //Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);

      //Set the solve flag
      Solve_which_system = Full_augmented;
     }
   }


 ///////////////////////////////////////////////////////////////////////
 // Non-inline functions for the BlockHopfLinearSolver class
 //////////////////////////////////////////////////////////////////////


 //======================================================================
 /// Clean up the memory that may have been allocated by the solver
 //=====================================================================
 BlockHopfLinearSolver::~BlockHopfLinearSolver()
 {
  if(A_pt!=0) {delete A_pt;}
  if(E_pt!=0) {delete E_pt;}
  if(G_pt!=0) {delete G_pt;}
 }

 //===================================================================
 /// Use a block factorisation to solve the augmented system
 /// associated with a Hopf bifurcation.
 //===================================================================
 void BlockHopfLinearSolver::solve(Problem* const &problem_pt,
                                   DoubleVector &result)
 {

  // if the result is setup then it should not be distributed
#ifdef PARANOID
  if (result.built())
   {
    if (result.distributed())
     {
      throw OomphLibError("The result vector must not be distributed",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
#endif

  //Locally cache the pointer to the handler.
  HopfHandler* handler_pt =
   static_cast<HopfHandler*>(problem_pt->assembly_handler_pt());

  //Find the number of dofs in the augmented problem
  unsigned n_dof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
  this->build_distribution(dist);

  //Firstly, let's calculate the derivative of the residuals wrt
  //the parameter
  DoubleVector dRdparam(this->distribution_pt(),0.0);

  const double FD_step = 1.0e-8;
  {
   //Cache pointer to the parameter (second last entry in the vector
   double *param_pt = &problem_pt->dof(n_dof-2);
   //Backup the parameter
   double old_var = *param_pt;;
   //Increment the parameter
   *param_pt += FD_step;
   problem_pt->actions_after_change_in_bifurcation_parameter();

   //Now calculate the new residuals
   problem_pt->get_residuals(dRdparam);

   //Now calculate the difference assume original residuals in resul
   for(unsigned n=0;n<n_dof;n++)
    {
     dRdparam[n] = (dRdparam[n] - result[n])/FD_step;
    }

   //Reset the parameter
   *param_pt = old_var;
   problem_pt->actions_after_change_in_bifurcation_parameter();
  }

  //Switch the handler to "standard" mode
  handler_pt->solve_standard_system();

  //Find out the number of dofs in the non-standard problem
  n_dof = problem_pt->ndof();

  // update the distribution pt
  dist.build(problem_pt->communicator_pt(),n_dof,false);
  this->build_distribution(dist);

  //Setup storage for temporary vector
    DoubleVector y1(this->distribution_pt(),0.0),
   alpha(this->distribution_pt(),0.0);

  //Allocate storage for A which can be used in the resolve
  if(A_pt!=0) {delete A_pt;}
  A_pt = new DoubleVector(this->distribution_pt(),0.0);

  //We are going to do resolves using the underlying linear solver
  Linear_solver_pt->enable_resolve();

  //Solve the first system Jy1 = R
  Linear_solver_pt->solve(problem_pt,y1);

  //This should have set the sign of the jacobian, store it
  int sign_of_jacobian = problem_pt->sign_of_jacobian();

  //Now let's get the appropriate bit of alpha
  for(unsigned n=0;n<n_dof;n++) {alpha[n] = dRdparam[n];}

  //Resolve to find A
  Linear_solver_pt->resolve(alpha,*A_pt);

  //Now set to the complex system
  handler_pt->solve_complex_system();

  // update the distribution
  dist.build(problem_pt->communicator_pt(),n_dof*2,false);
  this->build_distribution(dist);

  //Resize the stored vector G
  if(G_pt!=0) {delete G_pt;}
  G_pt = new DoubleVector(this->distribution_pt(),0.0);

  //Solve the first Zg = (Mz -My)
  Linear_solver_pt->solve(problem_pt,*G_pt);

  //This should have set the sign of the complex matrix's determinant, multiply
  sign_of_jacobian *= problem_pt->sign_of_jacobian();

  //We can now construct our multipliers
  //Prepare to scale
  double dof_length=0.0, a_length=0.0, y1_length=0.0;
  //Loop over the standard number of dofs
  for(unsigned n=0;n<n_dof;n++)
   {
    if(std::fabs(problem_pt->dof(n)) > dof_length)
     {dof_length = std::fabs(problem_pt->dof(n));}
    if(std::fabs((*A_pt)[n]) > a_length) {a_length = std::fabs((*A_pt)[n]);}
    if(std::fabs(y1[n]) > y1_length)
     {y1_length = std::fabs(y1[n]);}
   }

  double a_mult = dof_length/a_length;
  double y1_mult = dof_length/y1_length;
  a_mult += FD_step; y1_mult += FD_step;
  a_mult *= FD_step; y1_mult *= FD_step;

  //Local storage for the product terms
  Vector<double> Jprod_a(2*n_dof,0.0), Jprod_y1(2*n_dof,0.0);

  //Temporary storage
  Vector<double> rhs(2*n_dof);

  //find the number of elements
  unsigned long n_element = problem_pt->mesh_pt()->nelement();

  //Calculate the product of the jacobian matrices, etc
  for(unsigned long e=0;e<n_element;e++)
   {
    GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the ndofs in each element
    unsigned n_var = elem_pt->ndof();
    //Get the jacobian matrices
    DenseMatrix<double> jac(n_var), jac_a(n_var), jac_y1(n_var);
    DenseMatrix<double> M(n_var), M_a(n_var), M_y1(n_var);
    //Get unperturbed jacobian
    elem_pt->get_jacobian_and_mass_matrix(rhs,jac,M);

    //Backup the dofs
    Vector<double> dof_bac(n_var);
    //Perturb the original dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      dof_bac[n] = problem_pt->dof(eqn_number);
      //Pertub by vector A
      problem_pt->dof(eqn_number) += a_mult*(*A_pt)[eqn_number];
     }

    //Now get the new jacobian and mass matrix
    elem_pt->get_jacobian_and_mass_matrix(rhs,jac_a,M_a);

    //Perturb the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      problem_pt->dof(eqn_number) = dof_bac[n];
      //Pertub by vector y1
      problem_pt->dof(eqn_number) += y1_mult*y1[eqn_number];
     }

    //Now get the new jacobian and mass matrix
    elem_pt->get_jacobian_and_mass_matrix(rhs,jac_y1,M_y1);

    //Reset the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      problem_pt->dof(eqn_number) = dof_bac[n];
     }

    //OK, now work out the products
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      double prod_a1=0.0, prod_y11=0.0;
      double prod_a2=0.0, prod_y12=0.0;
      for(unsigned m=0;m<n_var;m++)
       {
        const unsigned unknown = elem_pt->eqn_number(m);
        const double y = handler_pt->Phi[unknown];
        const double z = handler_pt->Psi[unknown];
        const double omega = handler_pt->Omega;
        //Real part (first line)
        prod_a1 += (jac_a(n,m) - jac(n,m))*y
         + omega*(M_a(n,m) - M(n,m))*z;
        prod_y11 += (jac_y1(n,m) - jac(n,m))*y
         + omega*(M_y1(n,m) - M(n,m))*z;
        //Imag par (second line)
        prod_a2 += (jac_a(n,m) - jac(n,m))*z
         - omega*(M_a(n,m) - M(n,m))*y;
        prod_y12 += (jac_y1(n,m) - jac(n,m))*z
         - omega*(M_y1(n,m) - M(n,m))*y;
       }
      Jprod_a[eqn_number] += prod_a1/a_mult;
      Jprod_y1[eqn_number] += prod_y11/y1_mult;
      Jprod_a[n_dof + eqn_number] += prod_a2/a_mult;
      Jprod_y1[n_dof + eqn_number] += prod_y12/y1_mult;
     }
   }

   //The assumption here is still that the result has been set to the
   //residuals.
   for(unsigned n=0;n<2*n_dof;n++)
    {
     rhs[n] = result[n_dof+n] - Jprod_y1[n];
    }

   //Temporary storage
   DoubleVector y2(this->distribution_pt(),0.0);

   // DoubleVector for rhs
   DoubleVector rhs2(this->distribution_pt(),0.0);
   for (unsigned i = 0; i < 2*n_dof; i++)
    {
     rhs2[i] = rhs[i];
    }

   //Solve it
   Linear_solver_pt->resolve(rhs2,y2);

   //Assemble the next RHS
   for(unsigned n=0;n<2*n_dof;n++)
    {
     rhs[n] = dRdparam[n_dof+n] - Jprod_a[n];
    }

   //Resive the storage
   if(E_pt!=0) {delete E_pt;}
   E_pt = new DoubleVector(this->distribution_pt(),0.0);

   //Solve for the next RHS
   for (unsigned i = 0; i < 2*n_dof; i++)
    {
     rhs2[i] = rhs[i];
    }
   Linear_solver_pt->resolve(rhs2,*E_pt);

   //We can now calculate the final corrections
   //We need to work out a large number of dot products
   double dot_c = 0.0, dot_d = 0.0, dot_e = 0.0, dot_f = 0.0, dot_g=0.0;
   double dot_h = 0.0;

   for(unsigned n=0;n<n_dof;n++)
    {
     //Get the appopriate entry
     const double Cn = handler_pt->C[n];
     dot_c += Cn*y2[n];
     dot_d += Cn*y2[n_dof+n];
     dot_e += Cn*(*E_pt)[n];
     dot_f += Cn*(*E_pt)[n_dof+n];
     dot_g += Cn*(*G_pt)[n];
     dot_h += Cn*(*G_pt)[n_dof+n];
    }

   //Now we should be able to work out the corrections
   double denom = dot_e*dot_h - dot_g*dot_f;

   //Copy the previous residuals
   double R31 = result[3*n_dof], R32 = result[3*n_dof+1];
   //Delta parameter
   const double delta_param =
    ((R32 - dot_d)*dot_g - (R31 - dot_c)*dot_h)/denom;
   //Delta frequency
   const double delta_w =  -((R32 - dot_d) + dot_f*delta_param)/(dot_h);

   //Load into the result vector
   result[3*n_dof] = delta_param;
   result[3*n_dof+1] = delta_w;

   //The corrections to the null vector
   for(unsigned n=0;n<2*n_dof;n++)
    {
     result[n_dof + n] = y2[n] - (*E_pt)[n]*delta_param - (*G_pt)[n]*delta_w;
    }

   //Finally add the corrections to the unknowns
   for(unsigned n=0;n<n_dof;n++)
    {
     result[n] = y1[n] - (*A_pt)[n]*delta_param;
    }

   //The sign of the jacobian is the previous signs multiplied by the
   //sign of the denominator
   problem_pt->sign_of_jacobian() = sign_of_jacobian*
    static_cast<int>(std::fabs(denom)/denom);

   //Switch things to our full solver
   handler_pt->solve_full_system();

   //If we are not storing things, clear the results
   if(!Enable_resolve)
    {
     //We no longer need to store the matrix
     Linear_solver_pt->disable_resolve();
     delete A_pt; A_pt=0;
     delete E_pt; E_pt=0;
     delete G_pt; G_pt=0;
    }
   //Otherwise, also store the problem pointer
   else
    {
     Problem_pt = problem_pt;
    }
 }

 //=====================================================================
 /// Solve for two right hand sides, required for (efficient) continuation
 /// because otherwise we have to store the inverse of the jacobian
 /// and the complex equivalent ... nasty.
 //=====================================================================
 void BlockHopfLinearSolver::solve_for_two_rhs(Problem* const &problem_pt,
                                               DoubleVector &result,
                                               const DoubleVector &rhs2,
                                               DoubleVector &result2)
 {

  // if the result is setup then it should not be distributed
#ifdef PARANOID
  if (result.built())
   {
    if (result.distributed())
     {
      throw OomphLibError("The result vector must not be distributed",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
  if (result2.built())
   {
    if (result2.distributed())
     {
      throw OomphLibError("The result2 vector must not be distributed",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
#endif

  //Locally cache the pointer to the handler.
  HopfHandler* handler_pt =
   static_cast<HopfHandler*>(problem_pt->assembly_handler_pt());

  //Find the number of dofs in the augmented problem
  unsigned n_dof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
  this->build_distribution(dist);

  // if the result vector is not setup then rebuild with distribution = global
  if (!result.built())
   {
    result.build(this->distribution_pt(),0.0);
   }
  if (!result2.built())
   {
    result2.build(this->distribution_pt(),0.0);
   }

  //Firstly, let's calculate the derivative of the residuals wrt
  //the parameter
  DoubleVector dRdparam(this->distribution_pt(),0.0);

  const double FD_step = 1.0e-8;
  {
   //Cache pointer to the parameter (second last entry in the vector
   double *param_pt = &problem_pt->dof(n_dof-2);
   //Backup the parameter
   double old_var = *param_pt;;
   //Increment the parameter
   *param_pt += FD_step;
   problem_pt->actions_after_change_in_bifurcation_parameter();

   //Now calculate the new residuals
   problem_pt->get_residuals(dRdparam);

   //Now calculate the difference assume original residuals in resul
   for(unsigned n=0;n<n_dof;n++)
    {
     dRdparam[n] = (dRdparam[n] - result[n])/FD_step;
    }

   //Reset the parameter
   *param_pt = old_var;
   problem_pt->actions_after_change_in_bifurcation_parameter();
  }

  //Switch the handler to "standard" mode
  handler_pt->solve_standard_system();

  //Find out the number of dofs in the non-standard problem
  n_dof = problem_pt->ndof();

  //
  dist.build(problem_pt->communicator_pt(),n_dof,false);
  this->build_distribution(dist);

  //Setup storage for temporary vector
  DoubleVector y1(this->distribution_pt(),0.0),
   alpha(this->distribution_pt(),0.0),
   y1_resolve(this->distribution_pt(),0.0);

  //Allocate storage for A which can be used in the resolve
  if(A_pt!=0) {delete A_pt;}
  A_pt = new DoubleVector(this->distribution_pt(),0.0);

  //We are going to do resolves using the underlying linear solver
  Linear_solver_pt->enable_resolve();

  //Solve the first system Jy1 =
  Linear_solver_pt->solve(problem_pt,y1);

  //This should have set the sign of the jacobian, store it
  int sign_of_jacobian = problem_pt->sign_of_jacobian();

  //Now let's get the appropriate bit of alpha
  for(unsigned n=0;n<n_dof;n++) {alpha[n] = dRdparam[n];}

  //Resolve to find A
  Linear_solver_pt->resolve(alpha,*A_pt);

  //Get the solution for the second rhs
  for(unsigned n=0;n<n_dof;n++) {alpha[n] = rhs2[n];}

  //Resolve to find y1_resolve
  Linear_solver_pt->resolve(alpha,y1_resolve);

  //Now set to the complex system
  handler_pt->solve_complex_system();

  // rebuild the Distribution
  dist.build(problem_pt->communicator_pt(),n_dof*2,false);
  this->build_distribution(dist);

  //Resize the stored vector G
  if(G_pt!=0) {delete G_pt;}
  G_pt = new DoubleVector(this->distribution_pt(),0.0);

  //Solve the first Zg = (Mz -My)
  Linear_solver_pt->solve(problem_pt,*G_pt);

  //This should have set the sign of the complex matrix's determinant, multiply
  sign_of_jacobian *= problem_pt->sign_of_jacobian();

  //We can now construct our multipliers
  //Prepare to scale
  double dof_length=0.0, a_length=0.0, y1_length=0.0, y1_resolve_length=0.0;
  //Loop over the standard number of dofs
  for(unsigned n=0;n<n_dof;n++)
   {
    if(std::fabs(problem_pt->dof(n)) > dof_length)
     {dof_length = std::fabs(problem_pt->dof(n));}
    if(std::fabs((*A_pt)[n]) > a_length) {a_length = std::fabs((*A_pt)[n]);}
    if(std::fabs(y1[n]) > y1_length)
     {y1_length = std::fabs(y1[n]);}
    if(std::fabs(y1_resolve[n]) > y1_resolve_length)
     {y1_resolve_length = std::fabs(y1[n]);}
   }


  double a_mult = dof_length/a_length;
  double y1_mult = dof_length/y1_length;
  double y1_resolve_mult = dof_length/y1_resolve_length;
  a_mult += FD_step; y1_mult += FD_step; y1_resolve_mult += FD_step;
  a_mult *= FD_step; y1_mult *= FD_step; y1_resolve_mult *= FD_step;

  //Local storage for the product terms
  Vector<double> Jprod_a(2*n_dof,0.0), Jprod_y1(2*n_dof,0.0);
  Vector<double> Jprod_y1_resolve(2*n_dof,0.0);

  //Temporary storage
  Vector<double> rhs(2*n_dof);

  //find the number of elements
  unsigned long n_element = problem_pt->mesh_pt()->nelement();

  //Calculate the product of the jacobian matrices, etc
  for(unsigned long e = 0;e<n_element;e++)
   {
    GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the ndofs in each element
    unsigned n_var = elem_pt->ndof();
    //Get the jacobian matrices
    DenseMatrix<double> jac(n_var), jac_a(n_var), jac_y1(n_var),
     jac_y1_resolve(n_var);
    DenseMatrix<double> M(n_var), M_a(n_var), M_y1(n_var), M_y1_resolve(n_var);
    //Get unperturbed jacobian
    elem_pt->get_jacobian_and_mass_matrix(rhs,jac,M);

    //Backup the dofs
    Vector<double> dof_bac(n_var);
    //Perturb the original dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      dof_bac[n] = problem_pt->dof(eqn_number);
      //Pertub by vector A
      problem_pt->dof(eqn_number) += a_mult*(*A_pt)[eqn_number];
     }

    //Now get the new jacobian and mass matrix
    elem_pt->get_jacobian_and_mass_matrix(rhs,jac_a,M_a);

    //Perturb the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      problem_pt->dof(eqn_number) = dof_bac[n];
      //Pertub by vector y1
      problem_pt->dof(eqn_number) += y1_mult*y1[eqn_number];
     }

    //Now get the new jacobian and mass matrix
    elem_pt->get_jacobian_and_mass_matrix(rhs,jac_y1,M_y1);

     //Perturb the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      problem_pt->dof(eqn_number) = dof_bac[n];
      //Pertub by vector y1
      problem_pt->dof(eqn_number) += y1_resolve_mult*y1_resolve[eqn_number];
     }

    //Now get the new jacobian and mass matrix
    elem_pt->get_jacobian_and_mass_matrix(rhs,jac_y1_resolve,M_y1_resolve);

    //Reset the dofs
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      problem_pt->dof(eqn_number) = dof_bac[n];
     }

    //OK, now work out the products
    for(unsigned n=0;n<n_var;n++)
     {
      unsigned eqn_number = elem_pt->eqn_number(n);
      double prod_a1=0.0, prod_y11=0.0, prod_y1_resolve1 = 0.0;
      double prod_a2=0.0, prod_y12=0.0, prod_y1_resolve2 = 0.0;
      for(unsigned m=0;m<n_var;m++)
       {
        const unsigned unknown = elem_pt->eqn_number(m);
        const double y = handler_pt->Phi[unknown];
        const double z = handler_pt->Psi[unknown];
        const double omega = handler_pt->Omega;
        //Real part (first line)
        prod_a1 += (jac_a(n,m) - jac(n,m))*y
         + omega*(M_a(n,m) - M(n,m))*z;
        prod_y11 += (jac_y1(n,m) - jac(n,m))*y
         + omega*(M_y1(n,m) - M(n,m))*z;
        prod_y1_resolve1 += (jac_y1_resolve(n,m) - jac(n,m))*y
         + omega*(M_y1_resolve(n,m) - M(n,m))*z;
        //Imag par (second line)
        prod_a2 += (jac_a(n,m) - jac(n,m))*z
         - omega*(M_a(n,m) - M(n,m))*y;
        prod_y12 += (jac_y1(n,m) - jac(n,m))*z
         - omega*(M_y1(n,m) - M(n,m))*y;
        prod_y1_resolve2 += (jac_y1_resolve(n,m) - jac(n,m))*z
         - omega*(M_y1_resolve(n,m) - M(n,m))*y;
       }
      Jprod_a[eqn_number] += prod_a1/a_mult;
      Jprod_y1[eqn_number] += prod_y11/y1_mult;
      Jprod_y1_resolve[eqn_number] += prod_y1_resolve1/y1_resolve_mult;
      Jprod_a[n_dof + eqn_number] += prod_a2/a_mult;
      Jprod_y1[n_dof + eqn_number] += prod_y12/y1_mult;
      Jprod_y1_resolve[n_dof + eqn_number] += prod_y1_resolve2/y1_resolve_mult;
     }
   }

   //The assumption here is still that the result has been set to the
   //residuals.
   for(unsigned n=0;n<2*n_dof;n++)
    {
     rhs[n] = result[n_dof+n] - Jprod_y1[n];
    }

   //Temporary storage
   DoubleVector y2(this->distribution_pt(),0.0);

   //Solve it
   DoubleVector temp_rhs(this->distribution_pt(),0.0);
   for (unsigned i = 0; i < 2*n_dof; i++)
    {
     temp_rhs[i] = rhs[i];
    }
   Linear_solver_pt->resolve(temp_rhs,y2);

   //Assemble the next RHS
   for(unsigned n=0;n<2*n_dof;n++)
    {
     rhs[n] = dRdparam[n_dof+n] - Jprod_a[n];
    }

   //Resive the storage
   if(E_pt!=0) {delete E_pt;}
   E_pt = new DoubleVector(this->distribution_pt(),0.0);

   //Solve for the next RHS
   for (unsigned i = 0; i < 2* n_dof; i++)
    {
     temp_rhs[i] = rhs[i];
    }
   Linear_solver_pt->resolve(temp_rhs,*E_pt);

   //Assemble the next RHS
   for(unsigned n=0;n<2*n_dof;n++)
    {
     rhs[n] = rhs2[n_dof+n] - Jprod_y1_resolve[n];
    }

   DoubleVector y2_resolve(this->distribution_pt(),0.0);
   for (unsigned i = 0; i < 2* n_dof; i++)
    {
     temp_rhs[i] = rhs[i];
    }
   Linear_solver_pt->resolve(temp_rhs,y2_resolve);


   //We can now calculate the final corrections
   //We need to work out a large number of dot products
   double dot_c = 0.0, dot_d = 0.0, dot_e = 0.0, dot_f = 0.0, dot_g=0.0;
   double dot_h = 0.0;

   double dot_c_resolve = 0.0, dot_d_resolve = 0.0;

   for(unsigned n=0;n<n_dof;n++)
    {
     //Get the appopriate entry
     const double Cn = handler_pt->C[n];
     dot_c += Cn*y2[n];
     dot_d += Cn*y2[n_dof+n];
     dot_c_resolve += Cn*y2_resolve[n];
     dot_d_resolve += Cn*y2_resolve[n_dof+n];
     dot_e += Cn*(*E_pt)[n];
     dot_f += Cn*(*E_pt)[n_dof+n];
     dot_g += Cn*(*G_pt)[n];
     dot_h += Cn*(*G_pt)[n_dof+n];
    }

   //Now we should be able to work out the corrections
   double denom = dot_e*dot_h - dot_g*dot_f;

   //Copy the previous residuals
   double R31 = result[3*n_dof], R32 = result[3*n_dof+1];
   //Delta parameter
   const double delta_param =
    ((R32 - dot_d)*dot_g - (R31 - dot_c)*dot_h)/denom;
   //Delta frequency
   const double delta_w =  -((R32 - dot_d) + dot_f*delta_param)/(dot_h);

   //Corrections
   double R31_resolve = rhs2[3*n_dof], R32_resolve = rhs2[3*n_dof+1];
   //Delta parameter
   const double delta_param_resolve =
    ((R32_resolve - dot_d_resolve)*dot_g -
     (R31_resolve - dot_c_resolve)*dot_h)/denom;
   //Delta frequency
   const double delta_w_resolve =
    -((R32_resolve - dot_d_resolve) + dot_f*delta_param_resolve)/(dot_h);


   //Load into the result vector
   result[3*n_dof] = delta_param;
   result[3*n_dof+1] = delta_w;

   //The corrections to the null vector
   for(unsigned n=0;n<2*n_dof;n++)
    {
     result[n_dof + n] = y2[n] - (*E_pt)[n]*delta_param - (*G_pt)[n]*delta_w;
    }

   //Finally add the corrections to the unknowns
   for(unsigned n=0;n<n_dof;n++)
    {
     result[n] = y1[n] - (*A_pt)[n]*delta_param;
    }

    //Load into the result vector
   result2[3*n_dof] = delta_param_resolve;
   result2[3*n_dof+1] = delta_w_resolve;

   //The corrections to the null vector
   for(unsigned n=0;n<2*n_dof;n++)
    {
     result2[n_dof + n] = y2_resolve[n] - (*E_pt)[n]*delta_param_resolve
      - (*G_pt)[n]*delta_w_resolve;
    }

   //Finally add the corrections to the unknowns
   for(unsigned n=0;n<n_dof;n++)
    {
     result2[n] = y1_resolve[n] - (*A_pt)[n]*delta_param_resolve;
    }


   //The sign of the jacobian is the previous signs multiplied by the
   //sign of the denominator
   problem_pt->sign_of_jacobian() = sign_of_jacobian*
    static_cast<int>(std::fabs(denom)/denom);

   //Switch things to our full solver
   handler_pt->solve_full_system();

   //If we are not storing things, clear the results
   if(!Enable_resolve)
    {
     //We no longer need to store the matrix
     Linear_solver_pt->disable_resolve();
     delete A_pt; A_pt=0;
     delete E_pt; E_pt=0;
     delete G_pt; G_pt=0;
    }
   //Otherwise, also store the problem pointer
   else
    {
     Problem_pt = problem_pt;
    }

 }


 //======================================================================
 //Hack the re-solve to use the block factorisation
 //======================================================================
 void BlockHopfLinearSolver::resolve(const DoubleVector & rhs,
                                     DoubleVector &result)
 {
  throw OomphLibError("resolve() is not implemented for this solver",
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);

 }


 ///////////////////////////////////////////////////////////////////////
 // Non-inline functions for the HopfHandler class
 //////////////////////////////////////////////////////////////////////

 //====================================================================
 /// Constructor: Initialise the hopf handler,
 /// by setting initial guesses for Phi, Psi  and calculating Count.
 /// If the system changes, a new  handler must be constructed.
 //===================================================================
 HopfHandler::HopfHandler(Problem* const &problem_pt,
                          double* const &parameter_pt) :
  Solve_which_system(0), Parameter_pt(parameter_pt), Omega(0.0)
 {
  //Set the problem pointer
  Problem_pt = problem_pt;
  //Set the number of non-augmented degrees of freedom
  Ndof = problem_pt->ndof();

  // create the linear algebra distribution for this solver
  // currently only global (non-distributed) distributions are allowed
  LinearAlgebraDistribution* dist_pt = new
   LinearAlgebraDistribution(problem_pt->communicator_pt(),Ndof,false);

  //Resize the vectors of additional dofs
  Phi.resize(Ndof);
  Psi.resize(Ndof);
  C.resize(Ndof);
  Count.resize(Ndof,0);

  //Loop over all the elements in the problem
  unsigned n_element = problem_pt->mesh_pt()->nelement();
  for(unsigned e=0;e<n_element;e++)
   {
    GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the local freedoms in an element
    unsigned n_var = elem_pt->ndof();
    for(unsigned n=0;n<n_var;n++)
     {
      //Increase the associated global equation number counter
      ++Count[elem_pt->eqn_number(n)];
     }
   }

  //Calculate the value Phi by
  //solving the system JPhi = dF/dlambda

  //Locally cache the linear solver
  LinearSolver* const linear_solver_pt = problem_pt->linear_solver_pt();

  //Save the status before entry to this routine
  bool enable_resolve = linear_solver_pt->is_resolve_enabled();

  //We need to do a resolve
  linear_solver_pt->enable_resolve();

  //Storage for the solution
  DoubleVector x(dist_pt,0.0);

  //Solve the standard problem, we only want to make sure that
  //we factorise the matrix, if it has not been factorised. We shall
  //ignore the return value of x.
  linear_solver_pt->solve(problem_pt,x);

  //Get the vector dresiduals/dparameter
  problem_pt->get_derivative_wrt_global_parameter(parameter_pt,x);

  // Copy rhs vector into local storage so it doesn't get overwritten
  // if the linear solver decides to initialise the solution vector, say,
  // which it's quite entitled to do!
  DoubleVector input_x(x);

  //Now resolve the system with the new RHS and overwrite the solution
  linear_solver_pt->resolve(input_x,x);

  //Restore the storage status of the linear solver
  if (enable_resolve)
   {
    linear_solver_pt->enable_resolve();
   }
  else
   {
    linear_solver_pt->disable_resolve();
   }

  //Normalise the solution x
  double length = 0.0;
  for(unsigned n=0;n<Ndof;n++) {length += x[n]*x[n];}
  length = sqrt(length);

  //Now add the real part of the null space components to the problem
  //unknowns and initialise it all
  //This is dumb at the moment ... fix with eigensolver?
  for(unsigned n=0;n<Ndof;n++)
   {
    problem_pt->Dof_pt.push_back(&Phi[n]);
    C[n] = Phi[n] = -x[n]/length;
   }

  //Set the imaginary part so that the appropriate residual is
  //zero initially (eigensolvers)
  for(unsigned n=0;n<Ndof;n+=2)
   {
    //Make sure that we are not at the end of an array of odd length
    if(n!=Ndof-1)
     {
      Psi[n] = C[n+1];
      Psi[n+1] = -C[n];
     }
    //If it's odd set the final entry to zero
    else
     {
      Psi[n] = 0.0;
     }
   }

  //Next add the imaginary parts of the null space components to the problem
  for(unsigned n=0;n<Ndof;n++)
   {
    problem_pt->Dof_pt.push_back(&Psi[n]);
   }
  //Now add the parameter
  problem_pt->Dof_pt.push_back(parameter_pt);
  //Finally add the frequency
  problem_pt->Dof_pt.push_back(&Omega);

  // rebuild the dof dist
  Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                         Ndof*3+2,false);
  //Remove all previous sparse storage used during Jacobian assembly
  Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);

  // delete the dist_pt
  delete dist_pt;
 }

 //====================================================================
 /// Constructor: Initialise the hopf handler,
 /// by setting initial guesses for Phi, Psi, Omega  and calculating Count.
 /// If the system changes, a new  handler must be constructed.
 //===================================================================
 HopfHandler::HopfHandler(Problem* const &problem_pt,
                          double* const &parameter_pt,
                          const double &omega,
                          const DoubleVector &phi,
                          const DoubleVector &psi) :
  Solve_which_system(0), Parameter_pt(parameter_pt), Omega(omega)
 {
  //Set the problem pointer
  Problem_pt = problem_pt;
  //Set the number of non-augmented degrees of freedom
  Ndof = problem_pt->ndof();

  //Resize the vectors of additional dofs
  Phi.resize(Ndof);
  Psi.resize(Ndof);
  C.resize(Ndof);
  Count.resize(Ndof,0);

  //Loop over all the elements in the problem
  unsigned n_element = problem_pt->mesh_pt()->nelement();
  for(unsigned e=0;e<n_element;e++)
   {
    GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
    //Loop over the local freedoms in an element
    unsigned n_var = elem_pt->ndof();
    for(unsigned n=0;n<n_var;n++)
     {
      //Increase the associated global equation number counter
      ++Count[elem_pt->eqn_number(n)];
     }
   }

  //Normalise the guess for phi
  double length = 0.0;
  for(unsigned n=0;n<Ndof;n++) {length += phi[n]*phi[n];}
  length = sqrt(length);

  //Now add the real part of the null space components to the problem
  //unknowns and initialise it all
  for(unsigned n=0;n<Ndof;n++)
   {
    problem_pt->Dof_pt.push_back(&Phi[n]);
    C[n] = Phi[n] = phi[n]/length;
    Psi[n] = psi[n]/length;
   }

  //Next add the imaginary parts of the null space components to the problem
  for(unsigned n=0;n<Ndof;n++)
   {
    problem_pt->Dof_pt.push_back(&Psi[n]);
   }

  //Now add the parameter
  problem_pt->Dof_pt.push_back(parameter_pt);
  //Finally add the frequency
  problem_pt->Dof_pt.push_back(&Omega);

  // rebuild the Dof_distribution_pt
  Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                         Ndof*3+2,false);
  //Remove all previous sparse storage used during Jacobian assembly
  Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 }



 //=======================================================================
 /// Destructor return the problem to its original (non-augmented) state
 //=======================================================================
 HopfHandler::~HopfHandler()
 {
  //If we are using the block solver reset the problem's linear solver
  //to the original one
  BlockHopfLinearSolver* block_hopf_solver_pt =
   dynamic_cast<BlockHopfLinearSolver*>(Problem_pt->linear_solver_pt());
  if(block_hopf_solver_pt)
   {
    //Reset the problem's linear solver
    Problem_pt->linear_solver_pt() = block_hopf_solver_pt->linear_solver_pt();
    //Delete the block solver
    delete block_hopf_solver_pt;
   }
  //Now return the problem to its original size
  Problem_pt->Dof_pt.resize(Ndof);
  Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                         Ndof,false);
  //Remove all previous sparse storage used during Jacobian assembly
  Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 }


 //=============================================================
 ///Get the number of elemental degrees of freedom
 //=============================================================
 unsigned HopfHandler::ndof(GeneralisedElement* const &elem_pt)
 {
  unsigned raw_ndof = elem_pt->ndof();
  switch(Solve_which_system)
   {
    //Full augmented system
   case 0:
    return (3*raw_ndof + 2);
    break;
    //Standard non-augmented system
   case 1:
    return raw_ndof;
    break;
    //Complex system
   case 2:
    return (2*raw_ndof);
    break;

   default:
    throw OomphLibError("Solve_which_system can only be 0,1 or 2",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 }

 //=============================================================
 ///Get the global equation number of the local unknown
 //============================================================
 unsigned long HopfHandler::eqn_number(GeneralisedElement* const &elem_pt,
                                       const unsigned &ieqn_local)
 {
  //Get the raw value
  unsigned raw_ndof = elem_pt->ndof();
  unsigned long global_eqn;
  if(ieqn_local < raw_ndof)
   {
    global_eqn = elem_pt->eqn_number(ieqn_local);
   }
  else if(ieqn_local < 2*raw_ndof)
   {
    global_eqn = Ndof + elem_pt->eqn_number(ieqn_local - raw_ndof);
   }
  else if(ieqn_local < 3*raw_ndof)
   {
    global_eqn = 2*Ndof + elem_pt->eqn_number(ieqn_local - 2*raw_ndof);
   }
  else if(ieqn_local == 3*raw_ndof)
   {
    global_eqn = 3*Ndof;
   }
  else
   {
    global_eqn = 3*Ndof+1;
   }
  return global_eqn;
 }

 //==================================================================
 ///Get the residuals
 //=================================================================
 void HopfHandler::get_residuals(GeneralisedElement* const &elem_pt,
                                 Vector<double> &residuals)
 {
  //Should only call get residuals for the full system
  if(Solve_which_system==0)
   {
    //Need to get raw residuals and jacobian
    unsigned raw_ndof = elem_pt->ndof();

    DenseMatrix<double> jacobian(raw_ndof), M(raw_ndof);
    //Get the basic residuals, jacobian and mass matrix
    elem_pt->get_jacobian_and_mass_matrix(residuals,jacobian,M);

    //Initialise the pen-ultimate residual
    residuals[3*raw_ndof] = -1.0/
     (double)(Problem_pt->mesh_pt()->nelement());
    residuals[3*raw_ndof+1] = 0.0;

    //Now multiply to fill in the residuals
    for(unsigned i=0;i<raw_ndof;i++)
     {
      residuals[raw_ndof+i] = 0.0;
      residuals[2*raw_ndof+i] = 0.0;
      for(unsigned j=0;j<raw_ndof;j++)
       {
        unsigned global_unknown = elem_pt->eqn_number(j);
        //Real part
        residuals[raw_ndof+i] +=
         jacobian(i,j)*Phi[global_unknown] + Omega*M(i,j)*Psi[global_unknown];
        //Imaginary part
        residuals[2*raw_ndof+i] +=
         jacobian(i,j)*Psi[global_unknown] - Omega*M(i,j)*Phi[global_unknown];
       }
      //Get the global equation number
      unsigned global_eqn = elem_pt->eqn_number(i);

      //Real part
      residuals[3*raw_ndof] += (Phi[global_eqn]*C[global_eqn])/
       Count[global_eqn];
      //Imaginary part
      residuals[3*raw_ndof+1] += (Psi[global_eqn]*C[global_eqn])/
       Count[global_eqn];
     }
   }
  else
   {
    throw OomphLibError("Solve_which_system can only be 0",
OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);

   }
 }


 //===============================================================
 /// \short Calculate the elemental Jacobian matrix "d equation
 /// / d variable".
//==================================================================
 void HopfHandler::get_jacobian(GeneralisedElement* const &elem_pt,
                                Vector<double> &residuals,
                                DenseMatrix<double> &jacobian)
 {
  //The standard case
  if(Solve_which_system==0)
   {
    unsigned augmented_ndof = ndof(elem_pt);
    unsigned raw_ndof = elem_pt->ndof();

    //Get the basic residuals and jacobian
    DenseMatrix<double> M(raw_ndof);
    elem_pt->get_jacobian_and_mass_matrix(residuals,jacobian,M);
    //Now fill in the actual residuals
    get_residuals(elem_pt,residuals);

    //Now the jacobian appears in other entries
    for(unsigned n=0;n<raw_ndof;++n)
     {
      for(unsigned m=0;m<raw_ndof;++m)
       {
        jacobian(raw_ndof+n,raw_ndof+m) = jacobian(n,m);
        jacobian(raw_ndof+n,2*raw_ndof+m) = Omega*M(n,m);
        jacobian(2*raw_ndof+n,2*raw_ndof+m) = jacobian(n,m);
        jacobian(2*raw_ndof+n,raw_ndof+m) = -Omega*M(n,m);
        unsigned global_eqn = elem_pt->eqn_number(m);
        jacobian(raw_ndof+n,3*raw_ndof+1) += M(n,m)*Psi[global_eqn];
        jacobian(2*raw_ndof+n,3*raw_ndof+1) -= M(n,m)*Phi[global_eqn];
       }

      unsigned local_eqn = elem_pt->eqn_number(n);
      jacobian(3*raw_ndof,raw_ndof+n) = C[local_eqn]/Count[local_eqn];
      jacobian(3*raw_ndof+1,2*raw_ndof+n) = C[local_eqn]/Count[local_eqn];
     }

    const double FD_step = 1.0e-8;

    Vector<double> newres_p(augmented_ndof), newres_m(augmented_ndof);

    //Loop over the dofs
    for(unsigned n=0;n<raw_ndof;n++)
     {
      //Just do the x's
      unsigned long global_eqn = eqn_number(elem_pt,n);
      double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
      double init = *unknown_pt;
      *unknown_pt += FD_step;

      //Get the new residuals
      get_residuals(elem_pt,newres_p);

      //Reset
      *unknown_pt = init;

      //Subtract
      *unknown_pt -= FD_step;
      get_residuals(elem_pt,newres_m);

      for(unsigned m=0;m<raw_ndof;m++)
       {
        jacobian(raw_ndof+m,n) =
         (newres_p[raw_ndof+m] - residuals[raw_ndof+m])/(FD_step);
        jacobian(2*raw_ndof+m,n) =
         (newres_p[2*raw_ndof+m] - residuals[2*raw_ndof+m])/(FD_step);
       }
      //Reset the unknown
      *unknown_pt = init;
     }

    {
     //Now do the global parameter
     double* unknown_pt = Problem_pt->Dof_pt[3*Ndof];
     double init = *unknown_pt;
     *unknown_pt += FD_step;

     Problem_pt->actions_after_change_in_bifurcation_parameter();
     //Get the new residuals
     get_residuals(elem_pt,newres_p);

     //Reset
     *unknown_pt = init;

     //Subtract
     *unknown_pt -= FD_step;
     get_residuals(elem_pt,newres_m);

     for(unsigned m=0;m<augmented_ndof-2;m++)
      {
       jacobian(m,3*raw_ndof) =
        (newres_p[m] - residuals[m])/FD_step;
      }
     //Reset the unknown
     *unknown_pt = init;
     Problem_pt->actions_after_change_in_bifurcation_parameter();
    }
   } //End of standard case
  //Normal case
  else if(Solve_which_system==1)
   {
    //Just get the normal jacobian and residuals
    elem_pt->get_jacobian(residuals,jacobian);
   }
  //Otherwise the augmented complex case
  else if(Solve_which_system==2)
   {
    unsigned raw_ndof = elem_pt->ndof();

    //Get the basic residuals and jacobian
    DenseMatrix<double> M(raw_ndof);
    elem_pt->get_jacobian_and_mass_matrix(residuals,jacobian,M);

    //We now need to fill in the other blocks
    for(unsigned n=0;n<raw_ndof;n++)
     {
      for(unsigned m=0;m<raw_ndof;m++)
       {
        jacobian(n,raw_ndof+m) = Omega*M(n,m);
        jacobian(raw_ndof+n,m) = -Omega*M(n,m);
        jacobian(raw_ndof+n,raw_ndof+m) = jacobian(n,m);
       }
     }

    //Now overwrite to fill in the residuals
    //The decision take is to solve for the mass matrix multiplied
    //terms in the residuals because they require no additional
    //information to assemble.
    for(unsigned n=0;n<raw_ndof;n++)
     {
      residuals[n] = 0.0;
      residuals[raw_ndof+n] = 0.0;
      for(unsigned m=0;m<raw_ndof;m++)
       {
        unsigned global_unknown = elem_pt->eqn_number(m);
        //Real part
        residuals[n] += M(n,m)*Psi[global_unknown];
        //Imaginary part
        residuals[raw_ndof+n] -= M(n,m)*Phi[global_unknown];
       }
     }
   } //End of complex augmented case
  else
   {
    throw OomphLibError("Solve_which_system can only be 0,1 or 2",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 }



 //==================================================================
 ///Get the derivatives of the augmented residuals with respect to
 ///a parameter
 //=================================================================
 void HopfHandler::get_dresiduals_dparameter(
  GeneralisedElement* const &elem_pt,
  double* const &parameter_pt,  Vector<double> &dres_dparam)
 {
  //Should only call get residuals for the full system
  if(Solve_which_system==0)
   {
    //Need to get raw residuals and jacobian
    unsigned raw_ndof = elem_pt->ndof();

    DenseMatrix<double> djac_dparam(raw_ndof), dM_dparam(raw_ndof);
    //Get the basic residuals, jacobian and mass matrix
    elem_pt->get_djacobian_and_dmass_matrix_dparameter(
     parameter_pt,dres_dparam,djac_dparam,dM_dparam);

    //Initialise the pen-ultimate residual, which does not
    //depend on the parameter
    dres_dparam[3*raw_ndof] = 0.0;
    dres_dparam[3*raw_ndof+1] = 0.0;

   //Now multiply to fill in the residuals
    for(unsigned i=0;i<raw_ndof;i++)
     {
      dres_dparam[raw_ndof+i] = 0.0;
      dres_dparam[2*raw_ndof+i] = 0.0;
      for(unsigned j=0;j<raw_ndof;j++)
       {
        unsigned global_unknown = elem_pt->eqn_number(j);
        //Real part
        dres_dparam[raw_ndof+i] +=
         djac_dparam(i,j)*Phi[global_unknown] +
         Omega*dM_dparam(i,j)*Psi[global_unknown];
        //Imaginary part
        dres_dparam[2*raw_ndof+i] +=
         djac_dparam(i,j)*Psi[global_unknown] -
         Omega*dM_dparam(i,j)*Phi[global_unknown];
       }
     }
   }
  else
   {
    throw OomphLibError("Solve_which_system can only be 0",
OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);

   }
 }




 //========================================================================
 /// Overload the derivative of the residuals and jacobian
 /// with respect to a parameter so that it breaks because it should not
 /// be required
 //========================================================================
 void HopfHandler::get_djacobian_dparameter(
  GeneralisedElement* const &elem_pt,
  double* const &parameter_pt,
  Vector<double> &dres_dparam,
  DenseMatrix<double> &djac_dparam)
 {
  std::ostringstream error_stream;
  error_stream <<
   "This function has not been implemented because it is not required\n";
  error_stream << "in standard problems.\n";
  error_stream <<
   "If you find that you need it, you will have to implement it!\n\n";

  throw OomphLibError(error_stream.str(),
OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
 }


 //=====================================================================
 /// Overload the hessian vector product function so that
 /// it breaks because it should not be required
 //========================================================================
 void HopfHandler::get_hessian_vector_products(
  GeneralisedElement* const &elem_pt,
  Vector<double> const &Y,
  DenseMatrix<double> const &C,
  DenseMatrix<double> &product)
 {
  std::ostringstream error_stream;
  error_stream <<
   "This function has not been implemented because it is not required\n";
  error_stream << "in standard problems.\n";
  error_stream <<
   "If you find that you need it, you will have to implement it!\n\n";

  throw OomphLibError(error_stream.str(),
OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
 }


 //==========================================================================
 /// Return the eigenfunction(s) associated with the bifurcation that
 /// has been detected in bifurcation tracking problems
 //==========================================================================
 void HopfHandler::get_eigenfunction(
  Vector<DoubleVector> &eigenfunction)
 {
  //There is a real and imaginary part of the null vector
  eigenfunction.resize(2);
  LinearAlgebraDistribution dist(Problem_pt->communicator_pt(),Ndof,false);
  //Rebuild the vector
  eigenfunction[0].build(&dist,0.0);
  eigenfunction[1].build(&dist,0.0);
  //Set the value to be the null vector
  for(unsigned n=0;n<Ndof;n++)
   {
    eigenfunction[0][n] = Phi[n];
    eigenfunction[1][n] = Psi[n];
   }
 }


 //====================================================================
 /// Set to solve the standard (underlying jacobian)  system
 //===================================================================
 void HopfHandler::solve_standard_system()
 {
  if(Solve_which_system!=1)
   {
    Solve_which_system = 1;
    //Restrict the problem to the standard variables only
    Problem_pt->Dof_pt.resize(Ndof);
    Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                           Ndof,false);
    //Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
   }
 }

 //====================================================================
 /// Set to solve the complex (jacobian and mass matrix)  system
 //===================================================================
 void HopfHandler::solve_complex_system()
 {
  //If we were not solving the complex system resize the unknowns
  //accordingly
  if(Solve_which_system!=2)
   {
    Solve_which_system=2;
    //Resize to the first Ndofs (will work whichever system we were
    //solving before)
    Problem_pt->Dof_pt.resize(Ndof);
    //Add the first (real) part of the eigenfunction back into the problem
    for(unsigned n=0;n<Ndof;n++)
     {
      Problem_pt->Dof_pt.push_back(&Phi[n]);
     }
    Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                           Ndof*2,false);
    //Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
   }
 }


 //=================================================================
 /// Set to Solve full system system
 //=================================================================
 void HopfHandler::solve_full_system()
 {
  //If we are starting from another system
  if(Solve_which_system)
   {
    Solve_which_system = 0;

    //Resize to the first Ndofs (will work whichever system we were
    //solving before)
    Problem_pt->Dof_pt.resize(Ndof);
    //Add the additional unknowns back into the problem
    for(unsigned n=0;n<Ndof;n++)
     {
      Problem_pt->Dof_pt.push_back(&Phi[n]);
     }
    for(unsigned n=0;n<Ndof;n++)
     {
      Problem_pt->Dof_pt.push_back(&Psi[n]);
     }
    //Now add the parameter
    Problem_pt->Dof_pt.push_back(Parameter_pt);
    //Finally add the frequency
    Problem_pt->Dof_pt.push_back(&Omega);

    //
    Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                           3*Ndof+2,false);
    //Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
   }
 }
}