iterative_linear_solver.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// Copyright (C) 2006-2025 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//====================================================================
// The actual solve functions for Iterative solvers.
 
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// Oomph-lib includes
#include "iterative_linear_solver.h"
 
// Required to force_ get templated builds of iterative solvers for
// sumofmatrices class.
#include "sum_of_matrices.h"
 
 
namespace oomph
{
  //==================================================================
  /// Default preconditioner for iterative solvers: The base
  /// class for preconditioners is a fully functional (if trivial!)
  /// preconditioner.
  //==================================================================
  IdentityPreconditioner IterativeLinearSolver::Default_preconditioner;
 
 
  ///////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////
 
 
  //==================================================================
  /// Re-solve the system defined by the last assembled Jacobian
  /// and the rhs vector specified here. Solution is returned in
  /// the vector result.
  //==================================================================
  template<typename MATRIX>
  void BiCGStab<MATRIX>::resolve(const DoubleVector& rhs, DoubleVector& result)
  {
    // We are re-solving
    Resolving = true;
 
#ifdef PARANOID
    if (Matrix_pt == 0)
    {
      throw OomphLibError("No matrix was stored -- cannot re-solve",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Call linear algebra-style solver
    this->solve(Matrix_pt, rhs, result);
 
    // Reset re-solving flag
    Resolving = false;
  }
 
  //==================================================================
  /// Solver: Takes pointer to problem and returns the results vector
  /// which contains the solution of the linear system defined by
  /// the problem's fully assembled Jacobian and residual vector.
  //==================================================================
  template<typename MATRIX>
  void BiCGStab<MATRIX>::solve(Problem* const& problem_pt, DoubleVector& result)
  {
    // Initialise timer
#ifdef OOMPH_HAS_MPI
    double t_start = TimingHelpers::timer();
#else
    clock_t t_start = clock();
#endif
 
    // We're not re-solving
    Resolving = false;
 
    // Get rid of any previously stored data
    clean_up_memory();
 
    // Get Jacobian matrix in format specified by template parameter
    // and nonlinear residual vector
    Matrix_pt = new MATRIX;
    DoubleVector f;
    problem_pt->get_jacobian(f, *Matrix_pt);
 
 
    // We've made the matrix, we can delete it...
    Matrix_can_be_deleted = true;
 
    // Doc time for setup
#ifdef OOMPH_HAS_MPI
    double t_end = TimingHelpers::timer();
    Jacobian_setup_time = t_end - t_start;
#else
    clock_t t_end = clock();
    Jacobian_setup_time = double(t_end - t_start) / CLOCKS_PER_SEC;
#endif
 
    if (Doc_time)
    {
      oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                 << std::endl;
    }
 
    // set the distribution
    if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt))
    {
      // the solver has the same distribution as the matrix if possible
      this->build_distribution(
        dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt)
          ->distribution_pt());
    }
    else
    {
      // the solver has the same distribution as the RHS
      this->build_distribution(f.distribution_pt());
    }
 
    // Call linear algebra-style solver
    if ((result.built()) &&
        (!(*result.distribution_pt() == *this->distribution_pt())))
    {
      LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
      result.build(this->distribution_pt(), 0.0);
      this->solve_helper(Matrix_pt, f, result);
      result.redistribute(&temp_global_dist);
    }
    else
    {
      this->solve_helper(Matrix_pt, f, result);
    }
 
    // Kill matrix unless it's still required for resolve
    if (!Enable_resolve) clean_up_memory();
  };
 
 
  //==================================================================
  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
  /// and returns the solution of the linear system.
  /// Algorithm and variable names based on "Numerical Linear Algebra
  /// for High-Performance Computers" by Dongarra, Duff, Sorensen & van
  /// der Vorst. SIAM  (1998), page 185.
  //==================================================================
  template<typename MATRIX>
  void BiCGStab<MATRIX>::solve_helper(DoubleMatrixBase* const& matrix_pt,
                                      const DoubleVector& rhs,
                                      DoubleVector& solution)
  {
#ifdef PARANOID
    // check that the rhs vector is setup
    if (!rhs.built())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The vectors rhs must be setup";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // check that the matrix is square
    if (matrix_pt->nrow() != matrix_pt->ncol())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The matrix at matrix_pt must be square.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // check that the matrix and the rhs vector have the same nrow()
    if (matrix_pt->nrow() != rhs.nrow())
    {
      std::ostringstream error_message_stream;
      error_message_stream
        << "The matrix and the rhs vector must have the same number of rows.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // if the matrix is distributable then it too should have the same
    // communicator as the rhs vector
    DistributableLinearAlgebraObject* dist_matrix_pt =
      dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
    if (dist_matrix_pt != 0)
    {
      if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The matrix matrix_pt must have the same communicator as the "
             "vectors"
          << " rhs and result must have the same communicator";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // if the matrix is not distributable then it the rhs vector should not be
    // distributed
    else
    {
      if (rhs.distribution_pt()->distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The matrix (matrix_pt) is not distributable and therefore the rhs"
          << " vector must not be distributed";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // if the result vector is setup then check it has the same distribution
    // as the rhs
    if (solution.built())
    {
      if (!(*solution.distribution_pt() == *rhs.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream << "The solution vector distribution has been "
                                "setup; it must have the "
                             << "same distribution as the rhs vector.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // setup the solution if it is not
    if (!solution.distribution_pt()->built())
    {
      solution.build(this->distribution_pt(), 0.0);
    }
    // zero
    else
    {
      solution.initialise(0.0);
    }
 
    // Get number of dofs
    // unsigned n_dof=rhs.size();
    unsigned nrow_local = this->nrow_local();
 
    // Time solver
#ifdef OOMPH_HAS_MPI
    double t_start = TimingHelpers::timer();
#else
    clock_t t_start = clock();
#endif
 
    // Initialise: Zero initial guess so the initial residual is
    // equal to the RHS, i.e. the nonlinear residual
    DoubleVector residual(rhs);
    double residual_norm = residual.norm();
    double rhs_norm = residual_norm;
    if (rhs_norm == 0.0) rhs_norm = 1.0;
    DoubleVector x(rhs.distribution_pt(), 0.0);
 
    // Hat residual by copy operation
    DoubleVector r_hat(residual);
 
    // Normalised residual
    double normalised_residual_norm = residual_norm / rhs_norm;
 
    // if required will document convergence history to screen or file (if
    // stream open)
    if (Doc_convergence_history)
    {
      if (!Output_file_stream.is_open())
      {
        oomph_info << 0 << " " << normalised_residual_norm << std::endl;
      }
      else
      {
        Output_file_stream << 0 << " " << normalised_residual_norm << std::endl;
      }
    }
 
    // Check immediate convergence
    if (normalised_residual_norm < Tolerance)
    {
      if (Doc_time)
      {
        oomph_info << "BiCGStab converged immediately" << std::endl;
      }
      solution = x;
 
      // Doc time for solver
      double t_end = TimingHelpers::timer();
      Solution_time = t_end - t_start;
 
      if (Doc_time)
      {
        oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                   << std::endl;
      }
      return;
    }
 
    // Setup preconditioner only if we're not re-solving
    if (!Resolving)
    {
      // only setup the preconditioner if required
      if (Setup_preconditioner_before_solve)
      {
        // Setup preconditioner from the Jacobian matrix
        double t_start_prec = TimingHelpers::timer();
 
        preconditioner_pt()->setup(matrix_pt);
 
        // Doc time for setup of preconditioner
        double t_end_prec = TimingHelpers::timer();
        Preconditioner_setup_time = t_end_prec - t_start_prec;
 
        if (Doc_time)
        {
          oomph_info << "Time for setup of preconditioner  [sec]: "
                     << Preconditioner_setup_time << std::endl;
        }
      }
    }
    else
    {
      if (Doc_time)
      {
        oomph_info << "Setup of preconditioner is bypassed in resolve mode"
                   << std::endl;
      }
    }
 
    // Some auxiliary variables
    double rho_prev = 0.0;
    double alpha = 0.0;
    double omega = 0.0;
    double rho, beta, dot_prod, dot_prod_tt, dot_prod_ts;
    double s_norm, r_norm;
 
    // Some vectors
    DoubleVector p(this->distribution_pt(), 0.0),
      p_hat(this->distribution_pt(), 0.0), v(this->distribution_pt(), 0.0),
      z(this->distribution_pt(), 0.0), t(this->distribution_pt(), 0.0),
      s(this->distribution_pt(), 0.0);
 
    // Loop over max. number of iterations
    for (unsigned iter = 1; iter <= Max_iter; iter++)
    {
      // Dot product for rho
      rho = r_hat.dot(residual);
 
      // Breakdown?
      if (rho == 0.0)
      {
        oomph_info << "BiCGStab has broken down after " << iter << " iterations"
                   << std::endl;
        oomph_info << "Returning with current normalised residual of "
                   << normalised_residual_norm << std::endl;
      }
 
      // First step is different
      if (iter == 1)
      {
        p = residual;
      }
      else
      {
        beta = (rho / rho_prev) * (alpha / omega);
        for (unsigned i = 0; i < nrow_local; i++)
        {
          p[i] = residual[i] + beta * (p[i] - omega * v[i]);
        }
      }
 
      // Apply precondtitioner: p_hat=P^-1*p
      preconditioner_pt()->preconditioner_solve(p, p_hat);
 
      // Matrix vector product: v=A*p_hat
      matrix_pt->multiply(p_hat, v);
      dot_prod = r_hat.dot(v);
      alpha = rho / dot_prod;
      for (unsigned i = 0; i < nrow_local; i++)
      {
        s[i] = residual[i] - alpha * v[i];
      }
      s_norm = s.norm();
 
      // Normalised residual
      normalised_residual_norm = s_norm / rhs_norm;
 
      // if required will document convergence history to screen or file (if
      // stream open)
      if (Doc_convergence_history)
      {
        if (!Output_file_stream.is_open())
        {
          oomph_info << double(iter - 0.5) << " " << normalised_residual_norm
                     << std::endl;
        }
        else
        {
          Output_file_stream << double(iter - 0.5) << " "
                             << normalised_residual_norm << std::endl;
        }
      }
 
      // Converged?
      if (normalised_residual_norm < Tolerance)
      {
        for (unsigned i = 0; i < nrow_local; i++)
        {
          solution[i] = x[i] + alpha * p_hat[i];
        }
 
        if (Doc_time)
        {
          oomph_info << std::endl;
          oomph_info << "BiCGStab converged. Normalised residual norm: "
                     << normalised_residual_norm << std::endl;
          oomph_info << "Number of iterations to convergence: " << iter
                     << std::endl;
          oomph_info << std::endl;
        }
 
        // Store number of iterations taken
        Iterations = iter;
 
        // Doc time for solver
        double t_end = TimingHelpers::timer();
        Solution_time = t_end - t_start;
 
        if (Doc_time)
        {
          oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                     << std::endl;
        }
 
        return;
      }
 
      // Apply precondtitioner: z=P^-1*s
      preconditioner_pt()->preconditioner_solve(s, z);
 
      // Matrix vector product: t=A*z
      matrix_pt->multiply(z, t);
      dot_prod_ts = t.dot(s);
      dot_prod_tt = t.dot(t);
      omega = dot_prod_ts / dot_prod_tt;
      for (unsigned i = 0; i < nrow_local; i++)
      {
        x[i] += alpha * p_hat[i] + omega * z[i];
        residual[i] = s[i] - omega * t[i];
      }
      r_norm = residual.norm();
      rho_prev = rho;
 
      // Check convergence again
      normalised_residual_norm = r_norm / rhs_norm;
 
      // if required will document convergence history to screen or file (if
      // stream open)
      if (Doc_convergence_history)
      {
        if (!Output_file_stream.is_open())
        {
          oomph_info << iter << " " << normalised_residual_norm << std::endl;
        }
        else
        {
          Output_file_stream << iter << " " << normalised_residual_norm
                             << std::endl;
        }
      }
 
 
      if (normalised_residual_norm < Tolerance)
      {
        if (Doc_time)
        {
          oomph_info << std::endl;
          oomph_info << "BiCGStab converged. Normalised residual norm: "
                     << normalised_residual_norm << std::endl;
          oomph_info << "Number of iterations to convergence: " << iter
                     << std::endl;
          oomph_info << std::endl;
        }
        solution = x;
 
        // Store the number of itertions taken.
        Iterations = iter;
 
        // Doc time for solver
        double t_end = TimingHelpers::timer();
        Solution_time = t_end - t_start;
 
        if (Doc_time)
        {
          oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                     << std::endl;
        }
        return;
      }
 
 
      // Breakdown: Omega has to be >0 for to be able to continue
      if (omega == 0.0)
      {
        oomph_info << std::endl;
        oomph_info << "BiCGStab breakdown with omega=0.0. "
                   << "Normalised residual norm: " << normalised_residual_norm
                   << std::endl;
        oomph_info << "Number of iterations so far: " << iter << std::endl;
        oomph_info << std::endl;
        solution = x;
 
        // Store the number of itertions taken.
        Iterations = iter;
 
        // Doc time for solver
        double t_end = TimingHelpers::timer();
        Solution_time = t_end - t_start;
 
        if (Doc_time)
        {
          oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                     << std::endl;
        }
        return;
      }
 
 
    } // end of iteration loop
 
 
    // No convergence
    oomph_info << std::endl;
    oomph_info << "BiCGStab did not converge to required tolerance! "
               << std::endl;
    oomph_info << "Returning with normalised residual norm: "
               << normalised_residual_norm << std::endl;
    oomph_info << "after " << Max_iter << " iterations." << std::endl;
    oomph_info << std::endl;
 
    solution = x;
 
    // Doc time for solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
 
    if (Doc_time)
    {
      oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                 << std::endl;
    }
 
    if (Throw_error_after_max_iter)
    {
      std::string err = "Solver failed to converge and you requested an error";
      err += " on convergence failures.";
      throw OomphLibError(
        err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
  }
 
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
 
  //==================================================================
  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
  /// and returns the solution of the linear system.
  /// Algorithm and variable names based on "Matrix Computations,
  /// 2nd Ed." Golub & van Loan, John Hopkins University Press(1989),
  /// page 529.
  //==================================================================
  template<typename MATRIX>
  void CG<MATRIX>::solve_helper(DoubleMatrixBase* const& matrix_pt,
                                const DoubleVector& rhs,
                                DoubleVector& solution)
  {
#ifdef PARANOID
    // check that the rhs vector is setup
    if (!rhs.built())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The vectors rhs must be setup";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // check that the matrix is square
    if (matrix_pt->nrow() != matrix_pt->ncol())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The matrix at matrix_pt must be square.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // check that the matrix and the rhs vector have the same nrow()
    if (matrix_pt->nrow() != rhs.nrow())
    {
      std::ostringstream error_message_stream;
      error_message_stream
        << "The matrix and the rhs vector must have the same number of rows.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // if the matrix is distributable then it too should have the same
    // communicator as the rhs vector
    DistributableLinearAlgebraObject* dist_matrix_pt =
      dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
    if (dist_matrix_pt != 0)
    {
      if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The matrix matrix_pt must have the same communicator as the "
             "vectors"
          << " rhs and result must have the same communicator";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // if the matrix is not distributable then it the rhs vector should not be
    // distributed
    else
    {
      if (rhs.distribution_pt()->distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The matrix (matrix_pt) is not distributable and therefore the rhs"
          << " vector must not be distributed";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // if the result vector is setup then check it has the same distribution
    // as the rhs
    if (solution.built())
    {
      if (!(*solution.distribution_pt() == *rhs.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream << "The solution vector distribution has been "
                                "setup; it must have the "
                             << "same distribution as the rhs vector.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // setup the solution if it is not
    if (!solution.distribution_pt()->built())
    {
      solution.build(this->distribution_pt(), 0.0);
    }
    // zero
    else
    {
      solution.initialise(0.0);
    }
 
    // Get number of dofs
    // unsigned n_dof=rhs.size();
    unsigned nrow_local = this->nrow_local();
 
    // Initialise counter
    unsigned counter = 0;
 
    // Time solver
    double t_start = TimingHelpers::timer();
 
    // Initialise: Zero initial guess so the initial residual is
    // equal to the RHS
    DoubleVector x(this->distribution_pt(), 0.0);
    DoubleVector residual(rhs);
    double residual_norm = residual.norm();
    double rhs_norm = residual_norm;
    if (rhs_norm == 0.0) rhs_norm = 1.0;
 
    // Normalised residual
    double normalised_residual_norm = residual_norm / rhs_norm;
 
    // if required will document convergence history to screen or file (if
    // stream open)
    if (Doc_convergence_history)
    {
      if (!Output_file_stream.is_open())
      {
        oomph_info << 0 << " " << normalised_residual_norm << std::endl;
      }
      else
      {
        Output_file_stream << 0 << " " << normalised_residual_norm << std::endl;
      }
    }
 
    // Check immediate convergence
    if (normalised_residual_norm < Tolerance)
    {
      if (Doc_time)
      {
        oomph_info << "CG converged immediately" << std::endl;
      }
      solution = x;
 
      // Doc time for solver
      double t_end = TimingHelpers::timer();
      Solution_time = t_end - t_start;
 
      if (Doc_time)
      {
        oomph_info << "Time for solve with CG  [sec]: " << Solution_time
                   << std::endl;
      }
      return;
    }
 
 
    // Setup preconditioner only if we're not re-solving
    if (!Resolving)
    {
      // only setup the preconditioner if required
      if (Setup_preconditioner_before_solve)
      {
        // Setup preconditioner from the Jacobian matrix
        double t_start_prec = TimingHelpers::timer();
 
        preconditioner_pt()->setup(matrix_pt);
 
        // Doc time for setup of preconditioner
        double t_end_prec = TimingHelpers::timer();
        Preconditioner_setup_time = t_end_prec - t_start_prec;
 
        if (Doc_time)
        {
          oomph_info << "Time for setup of preconditioner  [sec]: "
                     << Preconditioner_setup_time << std::endl;
        }
      }
    }
    else
    {
      if (Doc_time)
      {
        oomph_info << "Setup of preconditioner is bypassed in resolve mode"
                   << std::endl;
      }
    }
 
 
    // Auxiliary vectors
    // Vector<double> z(n_dof),p(n_dof),jacobian_times_p(n_dof,0.0);
    DoubleVector z(this->distribution_pt(), 0.0),
      p(this->distribution_pt(), 0.0),
      jacobian_times_p(this->distribution_pt(), 0.0);
 
    // Auxiliary values
    double alpha, beta, rz;
    double prev_rz = 0.0;
 
    // Main iteration
    while ((normalised_residual_norm > Tolerance) && (counter != Max_iter))
    {
      // Apply precondtitioner: z=P^-1*r
      preconditioner_pt()->preconditioner_solve(residual, z);
 
      // P vector is computed differently for first and subsequent steps
      if (counter == 0)
      {
        p = z;
        rz = residual.dot(z);
      }
      // Subsequent steps
      else
      {
        rz = residual.dot(z);
        beta = rz / prev_rz;
        for (unsigned i = 0; i < nrow_local; i++)
        {
          p[i] = z[i] + beta * p[i];
        }
      }
 
 
      // Matrix vector product
      matrix_pt->multiply(p, jacobian_times_p);
      double pq = p.dot(jacobian_times_p);
      alpha = rz / pq;
 
      // Update
      prev_rz = rz;
      for (unsigned i = 0; i < nrow_local; i++)
      {
        x[i] += alpha * p[i];
        residual[i] -= alpha * jacobian_times_p[i];
      }
 
 
      // Calculate the 2norm
      residual_norm = residual.norm();
 
      // Difference between the initial and current 2norm residual
      normalised_residual_norm = residual_norm / rhs_norm;
 
 
      // if required will document convergence history to screen or file (if
      // stream open)
      if (Doc_convergence_history)
      {
        if (!Output_file_stream.is_open())
        {
          oomph_info << counter << " " << normalised_residual_norm << std::endl;
        }
        else
        {
          Output_file_stream << counter << " " << normalised_residual_norm
                             << std::endl;
        }
      }
 
      counter = counter + 1;
 
    } // end while
 
 
    if (counter >= Max_iter)
    {
      oomph_info << std::endl;
      oomph_info << "CG did not converge to required tolerance! " << std::endl;
      oomph_info << "Returning with normalised residual norm: "
                 << normalised_residual_norm << std::endl;
      oomph_info << "after " << counter << " iterations." << std::endl;
      oomph_info << std::endl;
    }
    else
    {
      if (Doc_time)
      {
        oomph_info << std::endl;
        oomph_info << "CG converged. Normalised residual norm: "
                   << normalised_residual_norm << std::endl;
        oomph_info << "Number of iterations to convergence: " << counter
                   << std::endl;
        oomph_info << std::endl;
      }
    }
 
 
    // Store number if iterations taken
    Iterations = counter;
 
    // Copy result back
    solution = x;
 
    // Doc time for solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
 
    if (Doc_time)
    {
      oomph_info << "Time for solve with CG  [sec]: " << Solution_time
                 << std::endl;
    }
 
    if ((counter >= Max_iter) && (Throw_error_after_max_iter))
    {
      std::string err = "Solver failed to converge and you requested an error";
      err += " on convergence failures.";
      throw OomphLibError(
        err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
 
  } // end CG
 
 
  //==================================================================
  /// Re-solve the system defined by the last assembled Jacobian
  /// and the rhs vector specified here. Solution is returned in
  /// the vector result.
  //==================================================================
  template<typename MATRIX>
  void CG<MATRIX>::resolve(const DoubleVector& rhs, DoubleVector& result)
  {
    // We are re-solving
    Resolving = true;
 
#ifdef PARANOID
    if (Matrix_pt == 0)
    {
      throw OomphLibError("No matrix was stored -- cannot re-solve",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Call linear algebra-style solver
    this->solve(Matrix_pt, rhs, result);
 
    // Reset re-solving flag
    Resolving = false;
  }
 
 
  //==================================================================
  /// Solver: Takes pointer to problem and returns the results vector
  /// which contains the solution of the linear system defined by
  /// the problem's fully assembled Jacobian and residual vector.
  //==================================================================
  template<typename MATRIX>
  void CG<MATRIX>::solve(Problem* const& problem_pt, DoubleVector& result)
  {
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // We're not re-solving
    Resolving = false;
 
    // Get rid of any previously stored data
    clean_up_memory();
 
    // Get Jacobian matrix in format specified by template parameter
    // and nonlinear residual vector
    Matrix_pt = new MATRIX;
    DoubleVector f;
    problem_pt->get_jacobian(f, *Matrix_pt);
 
    // We've made the matrix, we can delete it...
    Matrix_can_be_deleted = true;
 
    // Doc time for setup
    double t_end = TimingHelpers::timer();
    Jacobian_setup_time = t_end - t_start;
 
    if (Doc_time)
    {
      oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                 << std::endl;
    }
 
    // set the distribution
    if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt))
    {
      // the solver has the same distribution as the matrix if possible
      this->build_distribution(
        dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt)
          ->distribution_pt());
    }
    else
    {
      // the solver has the same distribution as the RHS
      this->build_distribution(f.distribution_pt());
    }
 
    // if the result vector is not setup
    if (!result.distribution_pt()->built())
    {
      result.build(this->distribution_pt(), 0.0);
    }
 
    // Call linear algebra-style solver
    if (!(*result.distribution_pt() == *this->distribution_pt()))
    {
      LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
      result.build(this->distribution_pt(), 0.0);
      this->solve_helper(Matrix_pt, f, result);
      result.redistribute(&temp_global_dist);
    }
    else
    {
      this->solve_helper(Matrix_pt, f, result);
    }
 
    // Kill matrix unless it's still required for resolve
    if (!Enable_resolve) clean_up_memory();
  };
 
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
 
  //==================================================================
  /// Self-test to be called inside solve_helper to ensure
  /// that all inputs are consistent and everything that needs to
  /// be built, is.
  //==================================================================
  template<typename MATRIX>
  void Smoother::check_validity_of_solve_helper_inputs(MATRIX* const& matrix_pt,
                                                       const DoubleVector& rhs,
                                                       DoubleVector& solution,
                                                       const double& n_dof)
  {
    // Check that if the matrix is not distributable then it should not
    // be distributed
    if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
    {
      if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
            ->distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream << "The matrix must not be distributed.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    } // if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)!=0)
    // Check that this rhs distribution is setup
    if (!rhs.built())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector distribution must be setup.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // Check that the rhs has the right number of global rows
    if (rhs.nrow() != n_dof)
    {
      std::ostringstream error_message_stream;
      error_message_stream << "RHS does not have the same dimension as "
                           << "the linear system";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // Check that the rhs is not distributed
    if (rhs.distribution_pt()->distributed())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector must not be distributed.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // Check that if the result is setup it matches the distribution
    // of the rhs
    if (solution.built())
    {
      if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream << "If the result distribution is setup then it "
                             << "must be the same as the rhs distribution";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    } // if (solution.built())
  } // End of check_validity_of_solve_helper_inputs
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //==================================================================
  /// Solver: Takes pointer to problem and returns the results
  /// vector which contains the solution of the linear system defined
  /// by the problem's fully assembled Jacobian and residual vector.
  //==================================================================
  template<typename MATRIX>
  void GS<MATRIX>::solve(Problem* const& problem_pt, DoubleVector& result)
  {
    // Reset the Use_as_smoother_flag as the solver is not being used
    // as a smoother
    Use_as_smoother = false;
 
    // Find the # of degrees of freedom (variables)
    unsigned n_dof = problem_pt->ndof();
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // We're not re-solving
    Resolving = false;
 
    // Get rid of any previously stored data
    clean_up_memory();
 
    // Set up the distribution
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
 
    // Assign the distribution to the LinearSolver
    this->build_distribution(dist);
 
    // Allocate space for the Jacobian matrix in format specified
    // by template parameter
    Matrix_pt = new MATRIX;
 
    // Get the nonlinear residuals vector
    DoubleVector f;
 
    // Assign the Jacobian and the residuals vector
    problem_pt->get_jacobian(f, *Matrix_pt);
 
    // We've made the matrix, we can delete it...
    Matrix_can_be_deleted = true;
 
    // Doc time for setup
    double t_end = TimingHelpers::timer();
    Jacobian_setup_time = t_end - t_start;
 
    // If time documentation is enabled
    if (Doc_time)
    {
      oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                 << std::endl;
    }
 
    // Call linear algebra-style solver
    this->solve_helper(Matrix_pt, f, result);
 
    // Kill matrix unless it's still required for resolve
    if (!Enable_resolve) clean_up_memory();
  } // End of solve
 
  //==================================================================
  /// Linear-algebra-type solver: Takes pointer to a matrix and
  /// rhs vector and returns the solution of the linear system.
  //==================================================================
  template<typename MATRIX>
  void GS<MATRIX>::solve_helper(DoubleMatrixBase* const& matrix_pt,
                                const DoubleVector& rhs,
                                DoubleVector& solution)
  {
    // Get number of dofs
    unsigned n_dof = rhs.nrow();
 
#ifdef PARANOID
    // Upcast the matrix to the appropriate type
    MATRIX* tmp_matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
 
    // PARANOID Run the self-tests to check the inputs are correct
    this->check_validity_of_solve_helper_inputs<MATRIX>(
      tmp_matrix_pt, rhs, solution, n_dof);
 
    // We don't need the pointer any more but we do still need the matrix
    // so just make tmp_matrix_pt a null pointer
    tmp_matrix_pt = 0;
#endif
 
    // Set up the solution if it is not
    if (!solution.distribution_pt()->built())
    {
      // Build it!
      solution.build(this->distribution_pt(), 0.0);
    }
    // If the solution is already set up
    else
    {
      // If we're inside the multigrid solver then as we traverse up the
      // hierarchy we use the smoother on the updated approximate solution.
      // As such, we should ONLY be resetting all the values to zero if
      // we're NOT inside the multigrid solver
      if (!Use_as_smoother)
      {
        // Initialise the vector with all entries set to zero
        solution.initialise(0.0);
      }
    } // if (!solution.distribution_pt()->built())
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // Copy the solution vector into x
    DoubleVector x(solution);
 
    // Create a vector to hold the residual. This will only be built if
    // we're not inside the multigrid solver
    DoubleVector current_residual;
 
    // Variable to hold the current residual norm. Only used if we're
    // not inside the multigrid solver
    double norm_res = 0.0;
 
    // Variables to hold the initial residual norm. Only used if we're
    // not inside the multigrid solver
    double norm_f = 0.0;
 
    // Initialise the value of Iterations
    Iterations = 0;
 
    // Calculate the residual only if we're not inside the multigrid solver
    if (!Use_as_smoother)
    {
      // Build the residual vector
      current_residual.build(this->distribution_pt(), 0.0);
 
      // Calculate the residual r=b-Ax
      matrix_pt->residual(x, rhs, current_residual);
 
      // Calculate the 2-norm of the residual vector
      norm_res = current_residual.norm();
 
      // Store the initial norm
      norm_f = norm_res;
 
      // If required, document the convergence history to screen or file (if
      // the output stream open)
      if (Doc_convergence_history)
      {
        // If the output file stream isn't open
        if (!Output_file_stream.is_open())
        {
          // Output the result to screen
          oomph_info << Iterations << " " << norm_res << std::endl;
        }
        // If the output file stream is open
        else
        {
          // Document the result to file
          Output_file_stream << Iterations << " " << norm_res << std::endl;
        }
      } // if (Doc_convergence_history)
    } // if (!Use_as_smoother)
 
    // Start of the main GS loop: run up to Max_iter times
    for (unsigned iter_num = 0; iter_num < Max_iter; iter_num++)
    {
      // Loop over rows
      for (unsigned i = 0; i < n_dof; i++)
      {
        double dummy = rhs[i];
        for (unsigned j = 0; j < i; j++)
        {
          dummy -= (*(matrix_pt))(i, j) * x[j];
        }
        for (unsigned j = (i + 1); j < n_dof; j++)
        {
          dummy -= (*(matrix_pt))(i, j) * x[j];
        }
        x[i] = dummy / (*(matrix_pt))(i, i);
      } // for(unsigned i=0;i<n_dof;i++)
 
      // Increment the value of Iterations
      Iterations++;
 
      // Calculate the residual only if we're not inside the multigrid solver
      if (!Use_as_smoother)
      {
        // Get residual
        matrix_pt->residual(x, rhs, current_residual);
 
        // Calculate the relative residual norm (i.e.
        // \frac{\|r_{i}\|}{\|r_{0}\|})
        norm_res = current_residual.norm() / norm_f;
 
        // If required will document convergence history to screen or file (if
        // stream open)
        if (Doc_convergence_history)
        {
          if (!Output_file_stream.is_open())
          {
            oomph_info << Iterations << " " << norm_res << std::endl;
          }
          else
          {
            Output_file_stream << Iterations << " " << norm_res << std::endl;
          }
        } // if (Doc_convergence_history)
 
        // Check the tolerance only if the residual norm is being computed
        if (norm_res < Tolerance)
        {
          // Break out of the for-loop
          break;
        }
      } // if (!Use_as_smoother)
    } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
 
    // Calculate the residual only if we're not inside the multigrid solver
    if (!Use_as_smoother)
    {
      // If time documentation is enabled
      if (Doc_time)
      {
        oomph_info << "\nGS converged. Residual norm: " << norm_res
                   << "\nNumber of iterations to convergence: " << Iterations
                   << "\n"
                   << std::endl;
      }
    } // if (!Use_as_smoother)
 
    // Copy result into result
    solution = x;
 
    // Doc. time for solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
    if (Doc_time)
    {
      oomph_info << "Time for solve with GS [sec]: " << Solution_time
                 << std::endl;
    }
 
    // If the solver failed to converge and the user asked for an error if
    // this happened
    if ((Iterations > Max_iter - 1) && (Throw_error_after_max_iter))
    {
      std::string error_message =
        "Solver failed to converge and you requested ";
      error_message += "an error on convergence failures.";
      throw OomphLibError(
        error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
  } // End of solve_helper
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //==================================================================
  /// Explicit template specialisation of the solver for CR
  /// matrices: Takes pointer to problem and returns the results
  /// vector which contains the solution of the linear system defined
  /// by the problem's fully assembled Jacobian and residual vector.
  //==================================================================
  void GS<CRDoubleMatrix>::solve(Problem* const& problem_pt,
                                 DoubleVector& result)
  {
    // Reset the Use_as_smoother_flag as the solver is not being used
    // as a smoother
    Use_as_smoother = false;
 
    // Find # of degrees of freedom (variables)
    unsigned n_dof = problem_pt->ndof();
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // We're not re-solving
    Resolving = false;
 
    // Get rid of any previously stored data
    clean_up_memory();
 
    // Set up the distribution
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
 
    // Build the distribution of the LinearSolver
    this->build_distribution(dist);
 
    // Allocate space for the Jacobian matrix in CRDoubleMatrix format
    Matrix_pt = new CRDoubleMatrix;
 
    // Build the matrix
    Matrix_pt->build(this->distribution_pt());
 
    // Allocate space for the residuals vector
    DoubleVector f;
 
    // Build it
    f.build(this->distribution_pt(), 0.0);
 
    // Get the global Jacobian and the accompanying residuals vector
    problem_pt->get_jacobian(f, *Matrix_pt);
 
    // Doc. time for setup
    double t_end = TimingHelpers::timer();
    Jacobian_setup_time = t_end - t_start;
 
    // Run the generic setup function
    setup_helper(Matrix_pt);
 
    // We've made the matrix, we can delete it...
    Matrix_can_be_deleted = true;
 
    // If time documentation is enabled
    if (Doc_time)
    {
      // Output the time for the assembly of the Jacobian to screen
      oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                 << std::endl;
    }
 
    // Call linear algebra-style solver
    solve_helper(Matrix_pt, f, result);
 
    // Kill matrix unless it's still required for resolve
    if (!Enable_resolve) clean_up_memory();
  } // End of solve
 
  //==================================================================
  /// Explicit template specialisation of the smoother_setup
  /// function for CR matrices. Set up the smoother for the matrix
  /// specified by the pointer. This definition of the smoother_setup
  /// has the added functionality that it sorts the entries in the
  /// given matrix so that the CRDoubleMatrix implementation of the
  /// solver can be used.
  //==================================================================
  void GS<CRDoubleMatrix>::setup_helper(DoubleMatrixBase* matrix_pt)
  {
    // Assume the matrix has been passed in from the outside so we must
    // not delete it. This is needed to avoid pre- and post-smoothers
    // deleting the same matrix in the MG solver. If this was originally
    // set to TRUE then this will be sorted out in the other functions
    // from which this was called
    Matrix_can_be_deleted = false;
 
    // Upcast the input matrix to system matrix to the type MATRIX
    Matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
    // The system matrix here is a CRDoubleMatrix. To make use of the
    // specific implementation of the solver for this type of matrix we
    // need to make sure the entries are arranged correctly
    Matrix_pt->sort_entries();
 
    // Now get access to the vector Index_of_diagonal_entries
    Index_of_diagonal_entries = Matrix_pt->get_index_of_diagonal_entries();
 
#ifdef PARANOID
    // Create a boolean variable to store the result of whether or not
    // the entries are in the correct order
    bool is_matrix_sorted = true;
 
    // Check to make sure the entries have been sorted properly
    is_matrix_sorted = Matrix_pt->entries_are_sorted();
 
    // If the entries are not sorted properly
    if (!is_matrix_sorted)
    {
      // Throw an error; the solver won't work unless the matrix is sorted
      // properly
      throw OomphLibError("Matrix is not sorted correctly",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
  } // End of setup_helper
 
  //==================================================================
  /// Explicit template specialisation of the solve_helper
  /// function for CR matrices. Exploiting the sparsity of the given
  /// matrix allows for a much faster iterative solver.
  //==================================================================
  void GS<CRDoubleMatrix>::solve_helper(DoubleMatrixBase* const& matrix_pt,
                                        const DoubleVector& rhs,
                                        DoubleVector& solution)
  {
    // Get number of dofs
    unsigned n_dof = rhs.nrow();
 
#ifdef PARANOID
    // Upcast the matrix to the appropriate type
    CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
    // PARANOID Run the self-tests to check the inputs are correct
    this->check_validity_of_solve_helper_inputs<CRDoubleMatrix>(
      cr_matrix_pt, rhs, solution, n_dof);
 
    // We don't need the pointer any more but we do still need the matrix
    // so just make tmp_matrix_pt a null pointer
    cr_matrix_pt = 0;
#endif
 
    // Setup the solution if it is not
    if (!solution.distribution_pt()->built())
    {
      solution.build(this->distribution_pt(), 0.0);
    }
    // If the solution is already set up
    else
    {
      // If we're inside the multigrid solver then as we traverse up the
      // hierarchy we use the smoother on the updated approximate solution.
      // As such, we should ONLY be resetting all the values to zero if
      // we're NOT inside the multigrid solver
      if (!Use_as_smoother)
      {
        // Initialise the vector with all entries set to zero
        solution.initialise(0.0);
      }
    } // if (!solution.distribution_pt()->built())
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // Copy the solution vector into x
    DoubleVector x(solution);
 
    // Create a vector to hold the residual. This will only be built if
    // we're not inside the multigrid solver
    DoubleVector current_residual;
 
    // Variable to hold the current residual norm. Only used if we're
    // not inside the multigrid solver
    double norm_res = 0.0;
 
    // Variables to hold the initial residual norm. Only used if we're
    // not inside the multigrid solver
    double norm_f = 0.0;
 
    // Initialise the value of Iterations
    Iterations = 0;
 
    // Calculate the residual only if we're not inside the multigrid solver
    if (!Use_as_smoother)
    {
      // Build the residual vector
      current_residual.build(this->distribution_pt(), 0.0);
 
      // Calculate the residual r=b-Ax
      matrix_pt->residual(x, rhs, current_residual);
 
      // Calculate the 2-norm of the residual vector
      norm_res = current_residual.norm();
 
      // Store the initial norm
      norm_f = norm_res;
 
      // If required will document convergence history to screen
      // or file (if stream is open)
      if (Doc_convergence_history)
      {
        if (!Output_file_stream.is_open())
        {
          oomph_info << Iterations << " " << norm_res << std::endl;
        }
        else
        {
          Output_file_stream << Iterations << " " << norm_res << std::endl;
        }
      } // if (!Use_as_smoother)
    } // if (Doc_convergence_history)
 
    // In each iteration of Gauss-Seidel we need to calculate
    //              x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
    // where L represents the strict lower triangular portion of A,
    // D represents the diagonal of A and U represents the strict upper
    // triangular portion of A. Since the value of inv(L+D)*rhs remains
    // constant throughout we shall calculate it now and store it to
    // avoid calculating it again with each iteration.
 
    // Create a temporary matrix pointer to allow for the extraction of the
    // row start, column index and value pointers
    CRDoubleMatrix* tmp_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
    // First acquire access to the value, row_start and column_index arrays
    // from the compressed row matrix
    const double* value_pt = tmp_matrix_pt->value();
    const int* row_start_pt = tmp_matrix_pt->row_start();
    const int* column_index_pt = tmp_matrix_pt->column_index();
 
    // We've finished using the temporary matrix pointer so make it a null
    // pointer
    tmp_matrix_pt = 0;
 
    // We can only calculate this constant term if the diagonal entries
    // of the matrix are nonzero so we check this first
    for (unsigned i = 0; i < n_dof; i++)
    {
      // Get the index of the last entry below or on the diagonal in row i
      unsigned diag_index = Index_of_diagonal_entries[i];
 
      if (unsigned(*(column_index_pt + diag_index)) != i)
      {
        std::string err_strng = "Gauss-Seidel cannot operate on a matrix with ";
        err_strng += "zero diagonal entries.";
        throw OomphLibError(
          err_strng, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    } // for (unsigned i=0;i<n_dof;i++)
 
    // If there were no zero diagonal entries we can calculate the vector
    //                 constant_term=inv(L+D)*rhs.
 
    // Create a vector to hold the constant term in the iteration
    DoubleVector constant_term(this->distribution_pt(), 0.0);
 
    // Copy the entries of rhs into the vector, constant_term
    constant_term = rhs;
 
    // Start by looping over the rows of constant_term
    for (unsigned i = 0; i < n_dof; i++)
    {
      // Get the index of the last entry below or on the diagonal in row i
      unsigned diag_index = Index_of_diagonal_entries[i];
 
      // Get the index of the first entry in row i
      unsigned row_i_start = *(row_start_pt + i);
 
      // If there are entries strictly below the diagonal then the
      // column index of the first entry in the i-th row cannot be
      // i (since we have ensured that nonzero entries exist on the
      // diagonal so this is the largest the column index of the
      // first entry in this row could be)
      if (unsigned(*(column_index_pt + row_i_start)) != i)
      {
        // Initialise the column index variable
        unsigned column_index = 0;
 
        // Loop over the entries below the diagonal on row i
        for (unsigned j = row_i_start; j < diag_index; j++)
        {
          // Find the column index of this nonzero entry
          column_index = *(column_index_pt + j);
 
          // Subtract the value of A(i,j)*constant_term(j) (where j<i)
          constant_term[i] -= (*(value_pt + j)) * constant_term[column_index];
        }
      } // if (*(column_index_pt+row_i_start)!=i)
 
      // Finish off by calculating constant_term(i)/A(i,i)
      constant_term[i] /= *(value_pt + diag_index);
    } // for (unsigned i=0;i<n_dof;i++)
 
    // Build a temporary vector to store the value of A*x with each iteration
    DoubleVector temp_vec(this->distribution_pt(), 0.0);
 
    // Outermost loop: Run max_iter times (the iteration number)
    for (unsigned iter_num = 0; iter_num < Max_iter; iter_num++)
    {
      // With each iteration we need to calculate
      //            x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
      //                 =-inv(L+D)*U*x_old+constant_term.
      // Since we've already calculated the vector constant_term, it remains
      // for us to calculate inv(L+D)*U*x. To do this we break the calculation
      // into two parts:
      //        (1) The matrix-vector multiplication temp_vec=U*x, and;
      //        (2) The matrix-vector multiplication inv(L+D)*temp_vec.
 
      //------------------------------
      // (1) Calculating temp_vec=U*x:
      //------------------------------
 
      // Auxiliary variable to store the index of the first entry on the
      // upper triangular portion of the matrix
      unsigned upper_tri_start = 0;
 
      // Auxiliary variable to store the index of the first entry in the
      // next row
      unsigned next_row_start = 0;
 
      // Set the temporary vector temp_vec to be initially be the zero vector
      temp_vec.initialise(0.0);
 
      // Loop over the rows of temp_vec (note: we do not loop over the final
      // row since all the entries on the last row of a strict upper triangular
      // matrix are zero)
      for (unsigned i = 0; i < n_dof - 1; i++)
      {
        // Get the index of the first entry on the upper triangular portion of
        // the matrix noting that the i-th entry in Index_of_diagonal_entries
        // will store the index of the diagonal entry of the i-th row
        upper_tri_start = Index_of_diagonal_entries[i] + 1;
 
        // Get the index of the first entry in the next row
        next_row_start = *(row_start_pt + i + 1);
 
        // Variable to store the column index of each entry
        unsigned column_index = 0;
 
        // Loop over all of the entries above the diagonal
        for (unsigned j = upper_tri_start; j < next_row_start; j++)
        {
          // Get the column index of this entry
          column_index = *(column_index_pt + j);
 
          // Update the value of temp_vec[i] by adding A(i,j)*x(j) (where j>i)
          temp_vec[i] += (*(value_pt + j)) * x[column_index];
        }
      } // for (unsigned i=0;i<n_dof;i++)
 
      //-----------------------------------
      // (2) Calculating inv(L+D)*temp_vec:
      //-----------------------------------
 
      // Now U*x=temp_vec has been calculated we need to calculate the vector
      //                 y=inv(L+D)*U*x=inv(L+D)*temp_vec.
      // Note: the i-th entry in y is given by
      //       y(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*y(j))/A(i,i).
      // We can see from this that to calculate the current entry in y we
      // use its previously calculated entries and the current entry in
      // temp_vec. As a result, we do not need two separate vectors y and
      // temp_vec for this calculation. Instead, we can just update the
      // entries in temp_vec (carefully!) using:
      //   temp_vec(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*temp_vec(j))/A(i,i).
 
      // Start by looping over the rows of rhs
      for (unsigned i = 0; i < n_dof; i++)
      {
        // Get the index of the last entry below or on the diagonal in row i
        unsigned diag_index = Index_of_diagonal_entries[i];
 
        // Get the index of the first entry in row i
        unsigned row_i_start = *(row_start_pt + i);
 
        // If there are no entries strictly below the diagonal then the
        // column index of the first entry in the i-th row will be greater
        // than or equal to i
        if (unsigned(*(column_index_pt + row_i_start)) != i)
        {
          // Initialise the column index variable
          unsigned column_index = 0;
 
          // Loop over the entries below the diagonal
          for (unsigned j = row_i_start; j < diag_index; j++)
          {
            // Find the column index of this nonzero entry
            column_index = *(column_index_pt + j);
 
            // Subtract the value of A(i,j)*y(j) (where j<i)
            temp_vec[i] -= (*(value_pt + j)) * temp_vec[column_index];
          }
        } // if (*(column_index_pt+row_i_start)!=i)
 
        // Finish off by calculating temp_vec(i)/A(i,i)
        temp_vec[i] /= *(value_pt + diag_index);
      } // for (unsigned i=0;i<n_dof;i++)
 
      //-------------------------
      // Updating the solution x:
      //-------------------------
 
      // The updated solution is given by
      //            x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
      //                 =-temp_vec+constant_term,
      // We have calculated both constant_term and temp_vec so we simply
      // need to calculate x now
 
      // Set x to be the vector, constant_term
      x = constant_term;
 
      // Update x to give the next iterate
      x -= temp_vec;
 
      // Increment the value of Iterations
      Iterations++;
 
      // Calculate the residual only if we're not inside the multigrid solver
      if (!Use_as_smoother)
      {
        // Get residual
        matrix_pt->residual(x, rhs, current_residual);
 
        // Calculate the relative residual norm (i.e.
        // \frac{\|r_{i}\|}{\|r_{0}\|})
        norm_res = current_residual.norm() / norm_f;
 
        // If required will document convergence history to screen or file (if
        // stream open)
        if (Doc_convergence_history)
        {
          if (!Output_file_stream.is_open())
          {
            oomph_info << iter_num << " " << norm_res << std::endl;
          }
          else
          {
            Output_file_stream << iter_num << " " << norm_res << std::endl;
          }
        } // if (Doc_convergence_history)
 
        // Check the tolerance only if the residual norm is being computed
        if (norm_res < Tolerance)
        {
          // Break out of the for-loop
          break;
        }
      } // if (!Use_as_smoother)
    } // for (unsigned iter_num=0;iter_num<max_iter;iter_num++)
 
    // Calculate the residual only if we're not inside the multigrid solver
    if (!Use_as_smoother)
    {
      // If time documentation is enabled
      if (Doc_time)
      {
        oomph_info << "\nGS converged. Residual norm: " << norm_res
                   << "\nNumber of iterations to convergence: " << Iterations
                   << "\n"
                   << std::endl;
      }
    } // if (!Use_as_smoother)
 
    // Copy the solution into the solution vector
    solution = x;
 
    // Document the time for the solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
 
    // If time documentation is enabled
    if (Doc_time)
    {
      oomph_info << "Time for solve with Gauss-Seidel [sec]: " << Solution_time
                 << std::endl;
    }
 
    // If the solver failed to converge and the user asked for an error if
    // this happened
    if ((Iterations > Max_iter - 1) && (Throw_error_after_max_iter))
    {
      std::string error_message =
        "Solver failed to converge and you requested ";
      error_message += "an error on convergence failures.";
      throw OomphLibError(
        error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
  } // End of solve_helper
 
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
 
  //==================================================================
  /// Solver: Takes pointer to problem and returns the results
  /// vector which contains the solution of the linear system defined
  /// by the problem's fully assembled Jacobian and residual vector.
  //==================================================================
  template<typename MATRIX>
  void DampedJacobi<MATRIX>::solve(Problem* const& problem_pt,
                                   DoubleVector& result)
  {
    // Reset the Use_as_smoother_flag as the solver is not being used
    // as a smoother
    Use_as_smoother = false;
 
    // Find the # of degrees of freedom (variables)
    unsigned n_dof = problem_pt->ndof();
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // We're not re-solving
    Resolving = false;
 
    // Get rid of any previously stored data
    clean_up_memory();
 
    // Set up the distribution
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
 
    // Assign the distribution to the LinearSolver
    this->build_distribution(dist);
 
    // Allocate space for the Jacobian matrix in format specified
    // by template parameter
    Matrix_pt = new MATRIX;
 
    // Get the nonlinear residuals vector
    DoubleVector f;
 
    // Assign the Jacobian and the residuals vector
    problem_pt->get_jacobian(f, *Matrix_pt);
 
    // Extract the diagonal entries of the matrix and store them
    extract_diagonal_entries(Matrix_pt);
 
    // We've made the matrix, we can delete it...
    Matrix_can_be_deleted = true;
 
    // Doc time for setup
    double t_end = TimingHelpers::timer();
    Jacobian_setup_time = t_end - t_start;
 
    // If time documentation is enabled
    if (Doc_time)
    {
      oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                 << std::endl;
    }
 
    // Call linear algebra-style solver
    solve_helper(Matrix_pt, f, result);
 
    // Kill matrix unless it's still required for resolve
    if (!Enable_resolve) clean_up_memory();
  } // End of solve
 
  //==================================================================
  /// Linear-algebra-type solver: Takes pointer to a matrix and
  /// rhs vector and returns the solution of the linear system.
  //==================================================================
  template<typename MATRIX>
  void DampedJacobi<MATRIX>::solve_helper(DoubleMatrixBase* const& matrix_pt,
                                          const DoubleVector& rhs,
                                          DoubleVector& solution)
  {
    // Get number of dofs
    unsigned n_dof = rhs.nrow();
 
#ifdef PARANOID
    // Upcast the matrix to the appropriate type
    MATRIX* tmp_matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
 
    // PARANOID Run the self-tests to check the inputs are correct
    this->check_validity_of_solve_helper_inputs<MATRIX>(
      tmp_matrix_pt, rhs, solution, n_dof);
 
    // We don't need the pointer any more but we do still need the matrix
    // so just make tmp_matrix_pt a null pointer
    tmp_matrix_pt = 0;
#endif
 
    // Setup the solution if it is not
    if (!solution.distribution_pt()->built())
    {
      // Build the distribution of the solution vector if it hasn't been done
      // yet
      solution.build(this->distribution_pt(), 0.0);
    }
    // If the solution has already been set up
    else
    {
      // If we're inside the multigrid solver then as we traverse up the
      // hierarchy we use the smoother on the updated approximate solution.
      // As such, we should ONLY be resetting all the values to zero if
      // we're NOT inside the multigrid solver
      if (!Use_as_smoother)
      {
        // Initialise the vector with all entries set to zero
        solution.initialise(0.0);
      }
    } // if (!solution.distribution_pt()->built())
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // Initial guess isn't necessarily zero (restricted solution from finer
    // grids) therefore x needs to be assigned values from the input.
    DoubleVector x(solution);
 
    // Create a vector to store the value of the constant term, omega*inv(D)*r
    DoubleVector constant_term(this->distribution_pt(), 0.0);
 
    // Calculate the constant term vector
    for (unsigned i = 0; i < n_dof; i++)
    {
      // Assign the i-th entry of constant_term
      constant_term[i] = Omega * Matrix_diagonal[i] * rhs[i];
    }
 
    // Create a vector to hold the residual. This will only be built if
    // we're not inside the multigrid solver
    DoubleVector local_residual;
 
    // Variable to store the 2-norm of the residual vector. Only used
    // if we are not working inside the MG solver
    double norm_res = 0.0;
 
    // Variables to hold the initial residual norm. Only used if we're
    // not inside the multigrid solver
    double norm_f = 0.0;
 
    // Initialise the value of Iterations
    Iterations = 0;
 
    // Calculate the residual only if we're not inside the multigrid solver
    if (!Use_as_smoother)
    {
      // Build the local residual vector
      local_residual.build(this->distribution_pt(), 0.0);
 
      // Calculate the residual vector
      matrix_pt->residual(x, rhs, local_residual);
 
      // Calculate the 2-norm
      norm_res = local_residual.norm();
 
      // Store the initial norm
      norm_f = norm_res;
 
      // If required will document convergence history to screen
      // or file (if stream is open)
      if (Doc_convergence_history)
      {
        if (!Output_file_stream.is_open())
        {
          oomph_info << Iterations << " " << norm_res << std::endl;
        }
        else
        {
          Output_file_stream << Iterations << " " << norm_res << std::endl;
        }
      } // if (Doc_convergence_history)
    } // if (!Use_as_smoother)
 
    // Create a temporary vector to store the value of A*e
    // on each iteration and another to store the residual
    // at the end of each iteration
    DoubleVector temp_vec(this->distribution_pt(), 0.0);
 
    // Outermost loop: Run up to Max_iter times (the iteration number)
    for (unsigned iter_num = 0; iter_num < Max_iter; iter_num++)
    {
      // Calculate A*e
      matrix_pt->multiply(x, temp_vec);
 
      // Loop over each degree of freedom and update
      // the current approximation
      for (unsigned idof = 0; idof < n_dof; idof++)
      {
        // Scale the idof'th entry of temp_vec
        // by omega/A(i,i)
        temp_vec[idof] *= Omega * Matrix_diagonal[idof];
      }
 
      // Update the value of e (in the system Ae=r)
      x += constant_term;
      x -= temp_vec;
 
      // Increment the value of Iterations
      Iterations++;
 
      // Calculate the residual only if we're not inside the multigrid solver
      if (!Use_as_smoother)
      {
        // Get residual
        matrix_pt->residual(x, rhs, local_residual);
 
        // Calculate the 2-norm of the residual r=b-Ax
        norm_res = local_residual.norm() / norm_f;
 
        // If required, this will document convergence history to
        // screen or file (if the stream is open)
        if (Doc_convergence_history)
        {
          if (!Output_file_stream.is_open())
          {
            oomph_info << Iterations << " " << norm_res << std::endl;
          }
          else
          {
            Output_file_stream << Iterations << " " << norm_res << std::endl;
          }
        } // if (Doc_convergence_history)
 
        // Check the tolerance only if the residual norm is being computed
        if (norm_res < Tolerance)
        {
          // Break out of the for-loop
          break;
        }
      } // if (!Use_as_smoother)
    } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
 
    // Calculate the residual only if we're not inside the multigrid solver
    if (!Use_as_smoother)
    {
      // If time documentation is enabled
      if (Doc_time)
      {
        oomph_info << "\nDamped Jacobi converged. Residual norm: " << norm_res
                   << "\nNumber of iterations to convergence: " << Iterations
                   << "\n"
                   << std::endl;
      }
    } // if (!Use_as_smoother)
 
    // Copy the solution into the solution vector
    solution = x;
 
    // Doc. time for solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
    if (Doc_time)
    {
      oomph_info << "Time for solve with damped Jacobi [sec]: " << Solution_time
                 << std::endl;
    }
 
    // If the solver failed to converge and the user asked for an error if
    // this happened
    if ((Iterations > Max_iter - 1) && (Throw_error_after_max_iter))
    {
      std::string error_message =
        "Solver failed to converge and you requested ";
      error_message += "an error on convergence failures.";
      throw OomphLibError(
        error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
  } // End of solve_helper
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //==================================================================
  /// Re-solve the system defined by the last assembled Jacobian
  /// and the rhs vector specified here. Solution is returned in
  /// the vector result.
  //==================================================================
  template<typename MATRIX>
  void GMRES<MATRIX>::resolve(const DoubleVector& rhs, DoubleVector& result)
  {
    // We are re-solving
    Resolving = true;
 
#ifdef PARANOID
    if (Matrix_pt == 0)
    {
      throw OomphLibError("No matrix was stored -- cannot re-solve",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Call linear algebra-style solver
    this->solve(Matrix_pt, rhs, result);
 
    // Reset re-solving flag
    Resolving = false;
  }
 
 
  //==================================================================
  /// Solver: Takes pointer to problem and returns the results vector
  /// which contains the solution of the linear system defined by
  /// the problem's fully assembled Jacobian and residual vector.
  //==================================================================
  template<typename MATRIX>
  void GMRES<MATRIX>::solve(Problem* const& problem_pt, DoubleVector& result)
  {
    // Find # of degrees of freedom (variables)
    unsigned n_dof = problem_pt->ndof();
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // We're not re-solving
    Resolving = false;
 
    // Get rid of any previously stored data
    clean_up_memory();
 
    // setup the distribution
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
    this->build_distribution(dist);
 
    // Get Jacobian matrix in format specified by template parameter
    // and nonlinear residual vector
    Matrix_pt = new MATRIX;
    DoubleVector f;
    if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt) != 0)
    {
      if (dynamic_cast<CRDoubleMatrix*>(Matrix_pt) != 0)
      {
        dynamic_cast<CRDoubleMatrix*>(Matrix_pt)->build(
          this->distribution_pt());
        f.build(this->distribution_pt(), 0.0);
      }
    }
    problem_pt->get_jacobian(f, *Matrix_pt);
 
    // We've made the matrix, we can delete it...
    Matrix_can_be_deleted = true;
 
    // Doc time for setup
    double t_end = TimingHelpers::timer();
    Jacobian_setup_time = t_end - t_start;
 
    if (Doc_time)
    {
      oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                 << std::endl;
    }
 
    // If we want to compute the gradient for the globally convergent
    // Newton method, then do it here
    if (Compute_gradient)
    {
      // Compute it
      Matrix_pt->multiply_transpose(f, Gradient_for_glob_conv_newton_solve);
      // Set the flag
      Gradient_has_been_computed = true;
    }
 
    // Call linear algebra-style solver
    // If the result distribution is wrong, then redistribute
    // before the solve and return to original distribution
    // afterwards
    if ((result.built()) &&
        (!(*result.distribution_pt() == *this->distribution_pt())))
    {
      LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
      result.build(this->distribution_pt(), 0.0);
      this->solve_helper(Matrix_pt, f, result);
      result.redistribute(&temp_global_dist);
    }
    // Otherwise just solve
    else
    {
      this->solve_helper(Matrix_pt, f, result);
    }
 
    // Kill matrix unless it's still required for resolve
    if (!Enable_resolve) clean_up_memory();
  };
 
 
  //==========================================================================
  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
  /// and returns the solution of the linear system.
  /// based on the algorithm presented in Templates for the
  /// Solution of Linear Systems: Building Blocks for Iterative
  /// Methods, Barrett,
  /// Berry et al, SIAM, 2006 and the implementation in the IML++ library :
  /// http://math.nist.gov/iml++/
  //==========================================================================
  template<typename MATRIX>
  void GMRES<MATRIX>::solve_helper(DoubleMatrixBase* const& matrix_pt,
                                   const DoubleVector& rhs,
                                   DoubleVector& solution)
  {
    // Get number of dofs
    unsigned n_dof = rhs.nrow();
 
#ifdef PARANOID
    // PARANOID check that if the matrix is distributable then it should not be
    // then it should not be distributed
    if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
    {
      if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
            ->distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream << "The matrix must not be distributed.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // PARANOID check that this rhs distribution is setup
    if (!rhs.built())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector distribution must be setup.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that the rhs has the right number of global rows
    if (rhs.nrow() != n_dof)
    {
      throw OomphLibError(
        "RHS does not have the same dimension as the linear system",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that the rhs is not distributed
    if (rhs.distribution_pt()->distributed())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector must not be distributed.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that if the result is setup it matches the distribution
    // of the rhs
    if (solution.built())
    {
      if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream << "If the result distribution is setup then it "
                                "must be the same as the "
                             << "rhs distribution";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Reset the time spent applying the preconditioner
    Preconditioner_application_time = 0.0;
 
    // Set up the solution if it is not
    if (!solution.built())
    {
      solution.build(this->distribution_pt(), 0.0);
    }
    // Otherwise initialise to zero
    else
    {
      solution.initialise(0.0);
    }
 
    // Time solver
    double t_start = TimingHelpers::timer();
 
    // Relative residual
    double resid;
 
    // iteration counter
    unsigned iter = 1;
 
    // if not using iteration restart set Restart to n_dof
    if (!Iteration_restart)
    {
      Restart = n_dof;
    }
 
    // initialise vectors
    Vector<double> s(Restart + 1, 0);
    Vector<double> cs(Restart + 1);
    Vector<double> sn(Restart + 1);
    DoubleVector w(this->distribution_pt(), 0.0);
 
    // Setup preconditioner only if we're not re-solving
    if (!Resolving)
    {
      // only setup the preconditioner before solve if require
      if (Setup_preconditioner_before_solve)
      {
        // Setup preconditioner from the Jacobian matrix
        double t_start_prec = TimingHelpers::timer();
 
        // do not setup
        preconditioner_pt()->setup(matrix_pt);
 
        // Doc time for setup of preconditioner
        double t_end_prec = TimingHelpers::timer();
        Preconditioner_setup_time = t_end_prec - t_start_prec;
 
        if (Doc_time)
        {
          oomph_info << "Time for setup of preconditioner  [sec]: "
                     << Preconditioner_setup_time << std::endl;
        }
      }
    }
    else
    {
      if (Doc_time)
      {
        oomph_info << "Setup of preconditioner is bypassed in resolve mode"
                   << std::endl;
      }
    }
 
    // solve b-Jx = Mr for r (assumes x = 0);
    DoubleVector r(this->distribution_pt(), 0.0);
    if (Preconditioner_LHS)
    {
      // Start the timer
      double t_start_prec = TimingHelpers::timer();
 
      // Apply the preconditioner
      preconditioner_pt()->preconditioner_solve(rhs, r);
 
      // Calculate the time taken for the preconditioner solve
      Preconditioner_application_time +=
        (TimingHelpers::timer() - t_start_prec);
    }
    else
    {
      r = rhs;
    }
    double normb = 0;
 
    // compute norm(r)
    double* r_pt = r.values_pt();
    for (unsigned i = 0; i < n_dof; i++)
    {
      normb += r_pt[i] * r_pt[i];
    }
    normb = sqrt(normb);
 
    // set beta (the initial residual)
    double beta = normb;
 
    // compute initial relative residual
    if (normb == 0.0) normb = 1;
    resid = beta / normb;
 
    // if required will document convergence history to screen or file (if
    // stream open)
    if (Doc_convergence_history)
    {
      if (!Output_file_stream.is_open())
      {
        oomph_info << 0 << " " << resid << std::endl;
      }
      else
      {
        Output_file_stream << 0 << " " << resid << std::endl;
      }
    }
 
    // if GMRES converges immediately
    if (resid <= Tolerance)
    {
      if (Doc_time)
      {
        oomph_info << "GMRES converged immediately. Normalised residual norm: "
                   << resid << std::endl;
      }
      // Doc time for solver
      double t_end = TimingHelpers::timer();
      Solution_time = t_end - t_start;
 
      if (Doc_time)
      {
        // Doc the time taken for the preconditioner applications
        oomph_info << "Time for all preconditioner applications [sec]: "
                   << Preconditioner_application_time
                   << "\n\nTime for solve with GMRES  [sec]: " << Solution_time
                   << std::endl;
      }
      return;
    }
 
 
    // initialise vector of orthogonal basis vectors (v) and upper hessenberg
    // matrix H
    // NOTE: for implementation purpose the upper hessenberg matrix indexes are
    // are swapped so the matrix is effectively transposed
    Vector<DoubleVector> v;
    v.resize(Restart + 1);
    Vector<Vector<double>> H(Restart + 1);
 
    // while...
    while (iter <= Max_iter)
    {
      // set zeroth basis vector v[0] to r/beta
      v[0].build(this->distribution_pt(), 0.0);
      double* v0_pt = v[0].values_pt();
      for (unsigned i = 0; i < n_dof; i++)
      {
        v0_pt[i] = r_pt[i] / beta;
      }
 
      //
      s[0] = beta;
 
      // inner iteration counter for restarted version
      unsigned iter_restart;
 
      // perform iterations
      for (iter_restart = 0; iter_restart < Restart && iter <= Max_iter;
           iter_restart++, iter++)
      {
        // resize next column of upper hessenberg matrix
        H[iter_restart].resize(iter_restart + 2);
 
        // solve Jv[i] = Mw for w
        {
          DoubleVector temp(this->distribution_pt(), 0.0);
          if (Preconditioner_LHS)
          {
            // solve Jv[i] = Mw for w
            matrix_pt->multiply(v[iter_restart], temp);
 
            // Start the timer
            double t_start_prec = TimingHelpers::timer();
 
            // Apply the preconditioner
            preconditioner_pt()->preconditioner_solve(temp, w);
 
            // Calculate the time taken for the preconditioner solve
            Preconditioner_application_time +=
              (TimingHelpers::timer() - t_start_prec);
          }
          else
          {
            // Start the timer
            double t_start_prec = TimingHelpers::timer();
 
            // w=JM^{-1}v by saad p270
            preconditioner_pt()->preconditioner_solve(v[iter_restart], temp);
 
            // Calculate the time taken for the preconditioner solve
            Preconditioner_application_time +=
              (TimingHelpers::timer() - t_start_prec);
 
            // Do a matrix multiplication
            matrix_pt->multiply(temp, w);
          }
        }
 
        //
        double* w_pt = w.values_pt();
        for (unsigned k = 0; k <= iter_restart; k++)
        {
          //
          H[iter_restart][k] = 0.0;
          double* vk_pt = v[k].values_pt();
          for (unsigned i = 0; i < n_dof; i++)
          {
            H[iter_restart][k] += w_pt[i] * vk_pt[i];
          }
 
          //
          for (unsigned i = 0; i < n_dof; i++)
          {
            w_pt[i] -= H[iter_restart][k] * vk_pt[i];
          }
        }
 
        //
        {
          double temp_norm_w = 0.0;
          for (unsigned i = 0; i < n_dof; i++)
          {
            temp_norm_w += w_pt[i] * w_pt[i];
          }
          temp_norm_w = sqrt(temp_norm_w);
          H[iter_restart][iter_restart + 1] = temp_norm_w;
        }
 
        //
        v[iter_restart + 1].build(this->distribution_pt(), 0.0);
        double* v_pt = v[iter_restart + 1].values_pt();
        for (unsigned i = 0; i < n_dof; i++)
        {
          v_pt[i] = w_pt[i] / H[iter_restart][iter_restart + 1];
        }
 
        //
        for (unsigned k = 0; k < iter_restart; k++)
        {
          apply_plane_rotation(
            H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
        }
        generate_plane_rotation(H[iter_restart][iter_restart],
                                H[iter_restart][iter_restart + 1],
                                cs[iter_restart],
                                sn[iter_restart]);
        apply_plane_rotation(H[iter_restart][iter_restart],
                             H[iter_restart][iter_restart + 1],
                             cs[iter_restart],
                             sn[iter_restart]);
        apply_plane_rotation(s[iter_restart],
                             s[iter_restart + 1],
                             cs[iter_restart],
                             sn[iter_restart]);
 
        // compute current residual
        beta = std::fabs(s[iter_restart + 1]);
 
        // compute relative residual
        resid = beta / normb;
 
        // if required will document convergence history to screen or file (if
        // stream open)
        if (Doc_convergence_history)
        {
          if (!Output_file_stream.is_open())
          {
            oomph_info << iter << " " << resid << std::endl;
          }
          else
          {
            Output_file_stream << iter << " " << resid << std::endl;
          }
        }
 
        // if required tolerance found
        if (resid < Tolerance)
        {
          // update result vector
          update(iter_restart, H, s, v, solution);
 
          // document convergence
          if (Doc_time)
          {
            oomph_info << std::endl;
            oomph_info << "GMRES converged (1). Normalised residual norm: "
                       << resid << std::endl;
            oomph_info << "Number of iterations to convergence: " << iter
                       << std::endl;
            oomph_info << std::endl;
          }
 
          // Doc time for solver
          double t_end = TimingHelpers::timer();
          Solution_time = t_end - t_start;
 
          Iterations = iter;
 
          if (Doc_time)
          {
            // Doc the time taken for the preconditioner applications
            oomph_info << "Time for all preconditioner applications [sec]: "
                       << Preconditioner_application_time
                       << "\n\nTime for solve with GMRES  [sec]: "
                       << Solution_time << std::endl;
          }
          return;
        }
      }
 
      // update
      if (iter_restart > 0) update((iter_restart - 1), H, s, v, solution);
 
      // solve Mr = (b-Jx) for r
      {
        DoubleVector temp(this->distribution_pt(), 0.0);
        matrix_pt->multiply(solution, temp);
        double* temp_pt = temp.values_pt();
        const double* rhs_pt = rhs.values_pt();
        for (unsigned i = 0; i < n_dof; i++)
        {
          temp_pt[i] = rhs_pt[i] - temp_pt[i];
        }
 
        if (Preconditioner_LHS)
        {
          // Start the timer
          double t_start_prec = TimingHelpers::timer();
 
          preconditioner_pt()->preconditioner_solve(temp, r);
 
          // Calculate the time taken for the preconditioner solve
          Preconditioner_application_time +=
            (TimingHelpers::timer() - t_start_prec);
        }
      }
 
      // compute current residual
      beta = 0.0;
      r_pt = r.values_pt();
      for (unsigned i = 0; i < n_dof; i++)
      {
        beta += r_pt[i] * r_pt[i];
      }
      beta = sqrt(beta);
 
      // if relative residual within tolerance
      resid = beta / normb;
      if (resid < Tolerance)
      {
        if (Doc_time)
        {
          oomph_info << std::endl;
          oomph_info << "GMRES converged (2). Normalised residual norm: "
                     << resid << std::endl;
          oomph_info << "Number of iterations to convergence: " << iter
                     << std::endl;
          oomph_info << std::endl;
        }
 
        // Doc time for solver
        double t_end = TimingHelpers::timer();
        Solution_time = t_end - t_start;
 
        Iterations = iter;
 
        if (Doc_time)
        {
          // Doc the time taken for the preconditioner applications
          oomph_info << "Time for all preconditioner applications [sec]: "
                     << Preconditioner_application_time
                     << "\n\nTime for solve with GMRES  [sec]: "
                     << Solution_time << std::endl;
        }
        return;
      }
    }
 
    // DRAIG: I want to delete this but GMRES doesn't compute the residuals
    // properly after this scope when it doesn't converge within the specified
    // number of iterations. Not sure why. If you're bored, try to figure
    // out why.
    {
      DoubleVector temp(this->distribution_pt(), 0.0);
      matrix_pt->multiply(solution, temp);
      temp *= -1.0;
      temp += rhs;
      resid = temp.norm() / normb;
    }
 
    // otherwise GMRES failed convergence
    oomph_info << std::endl;
    oomph_info << "GMRES did not converge to required tolerance! " << std::endl;
    oomph_info << "Returning with normalised residual norm: " << resid
               << std::endl;
    oomph_info << "after " << Max_iter << " iterations." << std::endl;
    oomph_info << std::endl;
 
 
    if (Throw_error_after_max_iter)
    {
      std::string err = "Solver failed to converge and you requested an error";
      err += " on convergence failures.";
      throw OomphLibError(
        err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
 
    return;
 
  } // End GMRES
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //==================================================================
  /// Solver: Takes pointer to problem and returns the results vector
  /// which contains the solution of the linear system defined by
  /// the problem's fully assembled Jacobian and residual vector.
  //==================================================================
  void AugmentedProblemGMRES::solve(Problem* const& problem_pt,
                                    DoubleVector& result)
  {
    // Find # of degrees of freedom (in the non-augmented problem)
    unsigned n_dof = problem_pt->ndof();
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // We're not re-solving
    Resolving = false;
 
    // Get rid of any previously stored data
    clean_up_memory();
 
    // Create a distribution
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
 
    // ...now build it
    this->build_distribution(dist);
 
    // Get the Jacobian matrix and the nonlinear residual vector
    Matrix_pt = new CRDoubleMatrix;
    DoubleVector f;
    if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt) != 0)
    {
      if (dynamic_cast<CRDoubleMatrix*>(Matrix_pt) != 0)
      {
        dynamic_cast<CRDoubleMatrix*>(Matrix_pt)->build(
          this->distribution_pt());
        f.build(this->distribution_pt(), 0.0);
      }
    }
    problem_pt->get_jacobian(f, *Matrix_pt);
 
    // We've made the matrix, we can delete it...
    Matrix_can_be_deleted = true;
 
    // Doc time for setup
    double t_end = TimingHelpers::timer();
    Jacobian_setup_time = t_end - t_start;
 
    if (Doc_time)
    {
      oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                 << std::endl;
    }
 
    // Reset the Schur complement scalar entry value
    Schur_complement_scalar = 1.0;
 
    // Call linear algebra-style solver. If the result distribution is wrong,
    // then redistribute before the solve and return to original distribution
    // afterwards
    if ((!(*result.distribution_pt() == *this->distribution_pt())) &&
        result.built())
    {
      // Create a temporary copy of the current result distribution
      LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
 
      // Re-build the result vector as an augmented vector
      result.build(dist, 0.0);
 
      // Call the auxilliary helper function to do the actual solve
      this->solve_helper(Matrix_pt, f, result);
 
      // Re-distribute result vector
      result.redistribute(&temp_global_dist);
    }
    // Otherwise just solve
    else
    {
      // Call the auxilliary helper function to do the actual solve
      this->solve_helper(Matrix_pt, f, result);
    }
 
    // Kill matrix unless it's still required for resolve
    if (!Enable_resolve) clean_up_memory();
  };
 
  //==================================================================
  /// Re-solve the system defined by the last assembled Jacobian
  /// and the rhs vector specified here. Solution is returned in
  /// the vector result.
  //==================================================================
  void AugmentedProblemGMRES::resolve(const DoubleVector& rhs,
                                      DoubleVector& result)
  {
    // We are re-solving
    Resolving = true;
 
#ifdef PARANOID
    if (Matrix_pt == 0)
    {
      throw OomphLibError("No matrix was stored -- cannot re-solve",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Call linear algebra-style solver
    this->solve(Matrix_pt, rhs, result);
 
    // Reset re-solving flag
    Resolving = false;
  }
 
 
  //=============================================================================
  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
  /// and returns the solution of the linear system.
  /// based on the algorithm presented in Templates for the
  /// Solution of Linear Systems: Building Blocks for Iterative Methods,
  /// Barrett, Berry et al, SIAM, 2006 and the implementation in the IML++
  /// library : http://math.nist.gov/iml++/
  //=============================================================================
  void AugmentedProblemGMRES::solve_helper(
    DoubleMatrixBase* const& input_matrix_pt,
    const DoubleVector& rhs,
    DoubleVector& lhs)
  {
    // Get number of dofs (in the non-augmented problem)
    unsigned n_dof = input_matrix_pt->nrow();
 
    // How many dofs in the augmented problem?
    unsigned n_aug_dof = n_dof + 1;
 
    // Get a handle to this Problems communicator
    OomphCommunicator* comm_pt = this->distribution_pt()->communicator_pt();
 
    // Set up the (non-augmented) distribution
    LinearAlgebraDistribution dist(comm_pt, n_dof, false);
 
    // Set up the (augmented) distribution
    LinearAlgebraDistribution aug_dist(comm_pt, n_aug_dof, false);
 
    // If the input vector doesn't have size N+1 then something has gone wrong
    if ((rhs.nrow() != n_dof) || (lhs.nrow() != n_dof))
    {
      // Create an output stream
      std::ostringstream error_message_stream;
 
      // Create the error message
      error_message_stream << "RHS vector has " << rhs.nrow()
                           << " rows and the "
                           << "LHS vector has " << rhs.nrow()
                           << " rows but\nthey "
                           << "should both have " << n_dof << " rows!"
                           << std::endl;
 
      // Throw an error
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
#ifdef PARANOID
    // PARANOID check that if the matrix is distributable then it should not be
    // then it should not be distributed
    if (dynamic_cast<DistributableLinearAlgebraObject*>(input_matrix_pt) != 0)
    {
      if (dynamic_cast<DistributableLinearAlgebraObject*>(input_matrix_pt)
            ->distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream << "The matrix must not be distributed.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // PARANOID check that this rhs distribution is setup
    if (!rhs.built())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector distribution must be setup.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that the rhs is not distributed
    if (rhs.distribution_pt()->distributed())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector must not be distributed.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that if the result is setup it matches the distribution
    // of the rhs
    if (lhs.built())
    {
      if (!(*rhs.distribution_pt() == *lhs.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream << "If the result distribution is setup then it "
                                "must be the same as the "
                             << "rhs distribution";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Set up the solution if it is not
    if (!lhs.built())
    {
      // Use the augmented distribution which has size N+1 (not the solvers
      // distribution which has size N)
      lhs.build(dist, 0.0);
    }
    // Otherwise initialise to zero
    else
    {
      lhs.initialise(0.0);
    }
 
    // Create a vector to store the augmented LHS vector entries
    DoubleVector solution(aug_dist, 0.0);
 
    // Time solver
    double t_start = TimingHelpers::timer();
 
    // Relative residual
    double resid;
 
    // Iteration counter
    unsigned iter = 1;
 
    // if not using iteration restart set Restart to n_aug_dof
    if (!Iteration_restart)
    {
      Restart = n_aug_dof;
    }
 
    // initialise vectors
    Vector<double> s(Restart + 1, 0);
    Vector<double> cs(Restart + 1);
    Vector<double> sn(Restart + 1);
    DoubleVector w(aug_dist, 0.0);
 
    // Set up the preconditioner only if we're not re-solving
    if (!Resolving)
    {
      // Only setup the preconditioner before solve if require
      if (Setup_preconditioner_before_solve)
      {
        // Setup preconditioner from the Jacobian matrix
        double t_start_prec = TimingHelpers::timer();
 
        // Do not setup
        preconditioner_pt()->setup(input_matrix_pt);
 
        // Doc time for setup of preconditioner
        double t_end_prec = TimingHelpers::timer();
        Preconditioner_setup_time = t_end_prec - t_start_prec;
 
        // If we're meant to document timings
        if (Doc_time)
        {
          oomph_info << "Time for setup of preconditioner [sec]: "
                     << Preconditioner_setup_time << std::endl;
        }
      } // if (Setup_preconditioner_before_solve)
    }
    else
    {
      // If we're meant to document timings
      if (Doc_time)
      {
        oomph_info << "Setup of preconditioner is bypassed in resolve mode"
                   << std::endl;
      }
    } // if (!Resolving)
 
    // The default system matrix is simply the one that was passed in
    CRDoubleMatrix* matrix_pt = dynamic_cast<CRDoubleMatrix*>(input_matrix_pt);
 
    // Storage for the time taken for all preconditioner applications
    double t_prec_application_total = 0.0;
 
    //---------------------------------------------------------------
    // Compute the scalar component of the diagonal Schur complement
    // preconditioner. NOTE: Only compute it if it is currently set
    // to the default value of unity. It's quite unlikely that it
    // will actually have value one (up to machine precision) in
    // practice. This check ensures that we only compute the value
    // during the solve function and not during the resolve.
    //---------------------------------------------------------------
    // If the Schur complement value has its default value of one
    if (std::fabs(Schur_complement_scalar - 1.0) < 1.0e-14)
    {
      // If the bordering vectors haven't actually been properly created
      if ((B_pt == 0) || (C_pt == 0))
      {
        // Throw an error
        throw OomphLibError("Bordering vectors have not been created!",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      // Make sure the bordering vectors have the right size
      else if ((B_pt->nrow() != n_dof) || (C_pt->nrow() != n_dof))
      {
        // Throw an error
        throw OomphLibError("Bordering vectors do not have the right size!",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Allocate storage for A^{-1}b
      DoubleVector a_inv_b(this->distribution_pt(), 0.0);
 
      // Compute A^{-1}b with the preconditioner approximation to A^{-1}
      preconditioner_pt()->preconditioner_solve(*B_pt, a_inv_b);
 
      // Now compute the scalar component of the Schur complement preconditioner
      Schur_complement_scalar = C_pt->dot(a_inv_b);
 
      // If the Schur complement scalar value is pratically zero or is NaN
      if ((std::fabs(Schur_complement_scalar) < 1.0e-14) ||
          std::isnan(Schur_complement_scalar))
      {
        // Reset it to have value 2.0 (not 1.0 otherwise we would end up
        // recomputing it every resolve) as we "should" be able to use
        // any scalar value and end up with the same number of distinct
        // eigenvalues for the preconditioned system).
        Schur_complement_scalar = 2.0;
      }
    } // if (std::fabs(Schur_complement_scalar-1.0)<1.0e-14)
 
    // Storage for the augmented RHS vector
    DoubleVector rhs_aug(aug_dist, 0.0);
 
    // Loop over the first n_dof entries of the RHS vector
    for (unsigned i = 0; i < n_dof; i++)
    {
      // Copy the i-th entry over
      rhs_aug[i] = rhs[i];
    }
 
    // Update the RHS entry associated with the augmented equation
    rhs_aug[n_dof] = (*Rhs_pt);
 
    // Storage for b-Jx=Mr for r (assumes x = 0);
    DoubleVector r(aug_dist, 0.0);
 
    // If we're using LHS preconditioning
    if (Preconditioner_LHS)
    {
      // Initialise timer
      double t_prec_application_start = TimingHelpers::timer();
 
      // Use the helper function to apply the Schur complement preconditioner
      apply_schur_complement_preconditioner(rhs_aug, r);
 
      // Doc time for setup of preconditioner
      double t_prec_application_end = TimingHelpers::timer();
 
      // Update the preconditioner application time total
      t_prec_application_total +=
        (t_prec_application_end - t_prec_application_start);
    }
    // If we're using RHS preconditioning
    else
    {
      // We don't need to do any preconditioning here; just copy the vector over
      r = rhs_aug;
    } // if (Preconditioner_LHS)
 
    // Compute the norm of the RHS vector
    double normb = r.norm();
 
    // Set beta (the initial residual)
    double beta = normb;
 
    // If the RHS vector is (numerically equivalent to) the zero vector
    if (normb == 0.0)
    {
      // Initialise the norm of b as unity
      normb = 1.0;
    }
 
    // Compute the relative residual
    resid = beta / normb;
 
    // If required will document convergence history to screen or file (if
    // stream open)
    if (Doc_convergence_history)
    {
      if (!Output_file_stream.is_open())
      {
        oomph_info << 0 << " " << resid << std::endl;
      }
      else
      {
        Output_file_stream << 0 << " " << resid << std::endl;
      }
    } // if (Doc_convergence_history)
 
    // If GMRES converges immediately
    if (resid <= Tolerance)
    {
      if (Doc_time)
      {
        oomph_info << "AugmentedProblemGMRES converged immediately. Normalised "
                   << "residual norm: " << resid << std::endl;
      }
 
      // Doc time for solver
      double t_end = TimingHelpers::timer();
      Solution_time = t_end - t_start;
 
      if (Doc_time)
      {
        oomph_info
          << "Time for all preconditioner applications [sec]: "
          << t_prec_application_total
          << "\nTotal time for solve with AugmentedProblemGMRES [sec]: "
          << Solution_time << std::endl;
      }
 
      // We're done so finish here
      return;
    }
 
    // initialise vector of orthogonal basis vectors (v) and upper hessenberg
    // matrix H
    // NOTE: for implementation purpose the upper hessenberg matrix indexes are
    // are swapped so the matrix is effectively transposed
    Vector<DoubleVector> v;
    v.resize(Restart + 1);
    Vector<Vector<double>> H(Restart + 1);
 
    // while...
    while (iter <= Max_iter)
    {
      // Copy the values of r into the zeroth basis vector v[0]
      v[0] = r;
 
      // Now scale it so the zeroth basis vector v[0] is r/beta
      v[0] /= beta;
 
      // Set the first entry of s
      s[0] = beta;
 
      // Counter for the inner iteration counter (restarted version)
      unsigned iter_restart = 0;
 
      // Start performing iterations
      for (iter_restart = 0; iter_restart < Restart && iter <= Max_iter;
           iter_restart++, iter++)
      {
        // Resize next column of upper hessenberg matrix
        H[iter_restart].resize(iter_restart + 2);
 
        // solve Jv[i] = Mw for w
        {
          // Allocate space for the product J*v[i]
          DoubleVector temp(aug_dist, 0.0);
 
          // If we're using LHS preconditioning
          if (Preconditioner_LHS)
          {
            // Compute J*v[i] where J is now the augmented system matrix
            augmented_matrix_multiply(matrix_pt, v[iter_restart], temp);
 
            // Initialise timer
            double t_prec_application_start = TimingHelpers::timer();
 
            // Use the helper function to apply the Schur complement
            // preconditioner
            apply_schur_complement_preconditioner(temp, w);
 
            // End timer
            double t_prec_application_end = TimingHelpers::timer();
 
            // Update the preconditioner application time total
            t_prec_application_total +=
              (t_prec_application_end - t_prec_application_start);
          }
          else
          {
            // Initialise timer
            double t_prec_application_start = TimingHelpers::timer();
 
            // Compute w=JM^{-1}v by saad p270
            apply_schur_complement_preconditioner(v[iter_restart], temp);
 
            // Doc time for setup of preconditioner
            double t_prec_application_end = TimingHelpers::timer();
 
            // Update the preconditioner application time total
            t_prec_application_total +=
              (t_prec_application_end - t_prec_application_start);
 
            // Calculate w=Jv_m where v_m=M^{-1}v
            augmented_matrix_multiply(matrix_pt, temp, w);
          }
        } // Solve Jv[i]=Mw for w
 
        // Get a pointer to the entries of w
        double* w_pt = w.values_pt();
 
        // Loop over the rows of the Hessenberg matrix
        for (unsigned k = 0; k <= iter_restart; k++)
        {
          // Compute the (k,iter_restart)-th entry of the Hessenberg matrix
          // (which is stored in its transposed form)
          H[iter_restart][k] = w.dot(v[k]);
 
          // Get a pointer to the entries of the k-th basis vector
          double* vk_pt = v[k].values_pt();
 
          // Loop over the entries of w
          for (unsigned i = 0; i < n_aug_dof; i++)
          {
            // Update the i-th entry of w
            w_pt[i] -= H[iter_restart][k] * vk_pt[i];
          }
        } // for (unsigned k=0;k<=iter_restart;k++)
 
        // Calculate the subdiagonal Hessenberg entry
        H[iter_restart][iter_restart + 1] = w.norm();
 
        // Copy the entries of w into the next basis vector
        v[iter_restart + 1] = w;
 
        // Now normalise the basis vector with the appropriate Hessenberg entry
        v[iter_restart + 1] /= H[iter_restart][iter_restart + 1];
 
        // Loop over the rows of the Hessenberg matrix
        for (unsigned k = 0; k < iter_restart; k++)
        {
          // Apply a Givens rotation to the appropriate entries of H
          apply_plane_rotation(
            H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
        }
 
        generate_plane_rotation(H[iter_restart][iter_restart],
                                H[iter_restart][iter_restart + 1],
                                cs[iter_restart],
                                sn[iter_restart]);
        apply_plane_rotation(H[iter_restart][iter_restart],
                             H[iter_restart][iter_restart + 1],
                             cs[iter_restart],
                             sn[iter_restart]);
        apply_plane_rotation(s[iter_restart],
                             s[iter_restart + 1],
                             cs[iter_restart],
                             sn[iter_restart]);
 
        // Compute the current residual
        beta = std::fabs(s[iter_restart + 1]);
 
        // compute relative residual
        resid = beta / normb;
 
        // DRAIG: Delete...
        // oomph_info << "Residual at iteration "
        //       << iter << ": " << resid << std::endl;
 
        // if required will document convergence history to screen or file (if
        // stream open)
        if (Doc_convergence_history)
        {
          if (!Output_file_stream.is_open())
          {
            oomph_info << iter << " " << resid << std::endl;
          }
          else
          {
            Output_file_stream << iter << " " << resid << std::endl;
          }
        } // if (Doc_convergence_history)
 
        // If required tolerance found
        if (resid < Tolerance)
        {
          // update result vector
          update(iter_restart, H, s, v, solution);
 
          // Loop over the first n_dof entries of the RHS vector
          for (unsigned i = 0; i < n_dof; i++)
          {
            // Copy the i-th entry over
            lhs[i] = solution[i];
          }
 
          // Update the bifurcation parameter
          (*X_pt) -= solution[n_dof];
 
          // If we're meant to document convergence
          if (Doc_time)
          {
            // Output the convergence to screen
            oomph_info << "\nAugmentedProblemGMRES converged (1). Normalised "
                       << "residual norm: " << resid
                       << "\nNumber of iterations "
                       << "to convergence: " << iter << "\n"
                       << std::endl;
          }
 
          // Doc time for solver
          double t_end = TimingHelpers::timer();
          Solution_time = t_end - t_start;
 
          Iterations = iter;
 
          // If we're meant to document convergence
          if (Doc_time)
          {
            oomph_info
              << "Time for all preconditioner applications [sec]: "
              << t_prec_application_total
              << "\nTotal time for solve with AugmentedProblemGMRES  [sec]: "
              << Solution_time << std::endl;
          }
 
          // We're done now so jump out here
          return;
        }
#ifdef PARANOID
        // If the tolerance wasn't met
        else
        {
          // Warn the user if the sub-diagonal Hessenberg entry was almost
          // zero...
          if (std::fabs(H[iter_restart][iter_restart + 1]) < 1.0e-12)
          {
            // Create an output stream
            std::ostringstream warning_message_stream;
 
            // Create the warning message
            warning_message_stream << "Hessenberg subdiagonal entry (="
                                   << H[iter_restart][iter_restart + 1]
                                   << ") almost zero at iteration "
                                   << iter_restart << "\nbut residual "
                                   << "is only " << resid << " which does "
                                   << "not meet the tolerance (" << Tolerance
                                   << ")!" << std::endl;
 
            // Warn the user
            OomphLibWarning(warning_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
          }
        } // if (resid<Tolerance)
#endif
      } // for (iter_restart=0;iter_restart<Restart&&iter<=Max_iter;...
 
      // If we're not on the first iteration
      if (iter_restart > 0)
      {
        // Do an update
        update((iter_restart - 1), H, s, v, solution);
      }
 
      // Solve Mr = (b-Jx) for r
      {
        DoubleVector temp(aug_dist, 0.0);
        augmented_matrix_multiply(matrix_pt, solution, temp);
        double* temp_pt = temp.values_pt();
        const double* rhs_pt = rhs.values_pt();
        for (unsigned i = 0; i < n_dof; i++)
        {
          // Update the i-th entry of temp to be b-Jx
          temp_pt[i] = rhs_pt[i] - temp_pt[i];
        }
 
        // If we're using LHS preconditioning
        if (Preconditioner_LHS)
        {
          // Initialise timer
          double t_prec_application_start = TimingHelpers::timer();
 
          // Use the helper function to apply the Schur complement
          // preconditioner
          apply_schur_complement_preconditioner(temp, r);
 
          // Doc time for setup of preconditioner
          double t_prec_application_end = TimingHelpers::timer();
 
          // Update the preconditioner application time total
          t_prec_application_total +=
            (t_prec_application_end - t_prec_application_start);
        }
      } // Solve Mr = (b-Jx) for r
 
      // Compute the current residual
      beta = r.norm();
 
      // if relative residual within tolerance
      resid = beta / normb;
      if (resid < Tolerance)
      {
        // If we're meant to document convergence
        if (Doc_time)
        {
          oomph_info << "\nAugmentedProblemGMRES converged (2). Normalised "
                     << "residual norm: " << resid << "\nNumber of iterations "
                     << "to convergence: " << iter << "\n"
                     << std::endl;
        }
 
        // Loop over the first n_dof entries of the RHS vector
        for (unsigned i = 0; i < n_dof; i++)
        {
          // Copy the i-th entry over
          lhs[i] = solution[i];
        }
 
        // Update the bifurcation parameter
        (*X_pt) -= solution[n_dof];
 
        // Doc time for solver
        double t_end = TimingHelpers::timer();
        Solution_time = t_end - t_start;
 
        Iterations = iter;
 
        // If we're meant to document convergence
        if (Doc_time)
        {
          oomph_info
            << "Time for all preconditioner applications [sec]: "
            << t_prec_application_total
            << "\nTotal time for solve with AugmentedProblemGMRES [sec]: "
            << Solution_time << std::endl;
        }
 
        return;
      } // if (resid<Tolerance)
    } // while (iter<=Max_iter)
 
    // DRAIG: Delete...
    {
      DoubleVector temp(aug_dist, 0.0);
      augmented_matrix_multiply(matrix_pt, solution, temp);
      temp *= -1.0;
      temp += rhs;
      resid = temp.norm() / normb;
    }
 
    // Loop over the first n_dof entries of the RHS vector
    for (unsigned i = 0; i < n_dof; i++)
    {
      // Copy the i-th entry over
      lhs[i] = solution[i];
    }
 
    // Update the bifurcation parameter
    (*X_pt) -= solution[n_dof];
 
    // Otherwise AugmentedProblemGMRES failed convergence
    oomph_info << "\nAugmentedProblemGMRES did not converge to required "
               << "tolerance!\nReturning with normalised residual norm: "
               << resid << "\nafter " << Max_iter << " iterations.\n"
               << std::endl;
 
    if (Throw_error_after_max_iter)
    {
      std::string error_message_string = "Solver failed to converge and you ";
      error_message_string += "requested an error on convergence failures.";
      throw OomphLibError(
        error_message_string, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
 
    return;
  } // End AugmentedProblemGMRES
 
 
  // Ensure build of required objects
  template class BiCGStab<CCDoubleMatrix>;
  template class BiCGStab<CRDoubleMatrix>;
  template class BiCGStab<DenseDoubleMatrix>;
 
  template class CG<CCDoubleMatrix>;
  template class CG<CRDoubleMatrix>;
  template class CG<DenseDoubleMatrix>;
 
  template class GS<CCDoubleMatrix>;
  template class GS<CRDoubleMatrix>;
  template class GS<DenseDoubleMatrix>;
 
  template class DampedJacobi<CCDoubleMatrix>;
  template class DampedJacobi<CRDoubleMatrix>;
  template class DampedJacobi<DenseDoubleMatrix>;
 
  template class GMRES<CCDoubleMatrix>;
  template class GMRES<CRDoubleMatrix>;
  template class GMRES<DenseDoubleMatrix>;
 
  // Solvers for SumOfMatrices class
  template class BiCGStab<SumOfMatrices>;
  template class CG<SumOfMatrices>;
  template class GS<SumOfMatrices>;
  template class GMRES<SumOfMatrices>;
} // namespace oomph