iterative_linear_solver.h

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//====================================================================
// This header defines a class for linear iterative solvers
 
// Include guards
#ifndef OOMPH_ITERATIVE_LINEAR_SOLVER_HEADER
#define OOMPH_ITERATIVE_LINEAR_SOLVER_HEADER
 
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
 
// oomph-lib headers
#include "matrices.h"
#include "problem.h"
#include "linear_solver.h"
#include "preconditioner.h"
 
 
namespace oomph
{
  //===========================================================================
  /// Base class for all linear iterative solvers.
  /// This merely defines standard interfaces for linear iterative solvers,
  /// so that different solvers can be used in a clean and transparent manner.
  //===========================================================================
  class IterativeLinearSolver : public LinearSolver
  {
  public:
    /// Constructor: Set (default) trivial preconditioner and set
    /// defaults for tolerance and max. number of iterations
    IterativeLinearSolver() : Preconditioner_setup_time(0)
    {
      // Set pointer to default preconditioner
      Preconditioner_pt = &Default_preconditioner;
 
      // Set default convergence tolerance
      Tolerance = 1.0e-8;
 
      // Set maximum number of iterations
      Max_iter = 100;
 
      // set default for document convergence history
      Doc_convergence_history = false;
 
      // set default
      Setup_preconditioner_before_solve = true;
 
      Throw_error_after_max_iter = false;
 
      // By default the iterative solver is not used as preconditioner
      Use_iterative_solver_as_preconditioner = false;
 
      // Indicates whether this is the first time we call the solve
      // method
      First_time_solve_when_used_as_preconditioner = true;
    }
 
    /// Broken copy constructor
    IterativeLinearSolver(const IterativeLinearSolver&) = delete;
 
    /// Broken assignment operator
    void operator=(const IterativeLinearSolver&) = delete;
 
    /// Destructor (empty)
    virtual ~IterativeLinearSolver() {}
 
    /// Access function to preconditioner
    Preconditioner*& preconditioner_pt()
    {
      return Preconditioner_pt;
    }
 
    /// Access function to preconditioner (const version)
    Preconditioner* const& preconditioner_pt() const
    {
      return Preconditioner_pt;
    }
 
    /// Access to convergence tolerance
    double& tolerance()
    {
      return Tolerance;
    }
 
    /// Access to max. number of iterations
    unsigned& max_iter()
    {
      return Max_iter;
    }
 
    /// Number of iterations taken
    virtual unsigned iterations() const = 0;
 
    /// Enable documentation of the convergence history
    void enable_doc_convergence_history()
    {
      Doc_convergence_history = true;
    }
 
    /// Disable documentation of the convergence history
    void disable_doc_convergence_history()
    {
      Doc_convergence_history = false;
    }
 
    /// Write convergence history into file with specified filename
    /// (automatically switches on doc). Optional second argument is a string
    /// that can be used (as a zone title) to identify what case
    /// we're running (e.g. what combination of linear solver and
    /// preconditioner or parameter values are used).
    void open_convergence_history_file_stream(
      const std::string& file_name, const std::string& zone_title = "")
    {
      // start docing
      Doc_convergence_history = true;
 
      // Close if it's open
      if (Output_file_stream.is_open())
      {
        Output_file_stream.close();
      }
 
      // Open new one
      Output_file_stream.open(file_name.c_str());
 
      // Write tecplot zone header
      Output_file_stream << "VARIABLES=\"iterations\", \"scaled residual\""
                         << std::endl;
      Output_file_stream << "ZONE T=\"" << zone_title << "\"" << std::endl;
      Output_file_stream << 0 << " " << 1.0 << std::endl;
    }
 
    /// Close convergence history output stream
    void close_convergence_history_file_stream()
    {
      if (Output_file_stream.is_open()) Output_file_stream.close();
    }
 
    ///  returns the time taken to assemble the jacobian matrix and
    /// residual vector
    double jacobian_setup_time() const
    {
      return Jacobian_setup_time;
    }
 
    /// return the time taken to solve the linear system
    double linear_solver_solution_time() const
    {
      return Solution_time;
    }
 
    /// returns the the time taken to setup the preconditioner
    virtual double preconditioner_setup_time() const
    {
      return Preconditioner_setup_time;
    }
 
    /// Setup the preconditioner before the solve
    void enable_setup_preconditioner_before_solve()
    {
      Setup_preconditioner_before_solve = true;
    }
 
    /// Don't set up the preconditioner before the solve
    void disable_setup_preconditioner_before_solve()
    {
      Setup_preconditioner_before_solve = false;
    }
 
    /// Throw an error if we don't converge within max_iter
    void enable_error_after_max_iter()
    {
      Throw_error_after_max_iter = true;
    }
 
    /// Don't throw an error if we don't converge within max_iter (default).
    void disable_error_after_max_iter()
    {
      Throw_error_after_max_iter = false;
    }
 
    /// Enables the iterative solver be used as preconditioner (when
    /// calling the solve method it bypass the setup solver method --
    /// currently only used by Trilinos solver ---)
    void enable_iterative_solver_as_preconditioner()
    {
      Use_iterative_solver_as_preconditioner = true;
    }
 
    /// Disables the iterative solver be used as preconditioner (when
    /// calling the solve method it bypass the setup solver method --
    /// currently only used by Trilinos solver ---)
    void disable_iterative_solver_as_preconditioner()
    {
      Use_iterative_solver_as_preconditioner = false;
    }
 
  protected:
    /// Flag indicating if the convergence history is to be
    /// documented
    bool Doc_convergence_history;
 
    /// Output file stream for convergence history
    std::ofstream Output_file_stream;
 
    /// Default preconditioner:  The base
    /// class for preconditioners is a fully functional (if trivial!)
    /// preconditioner.
    static IdentityPreconditioner Default_preconditioner;
 
    /// Convergence tolerance
    double Tolerance;
 
    /// Maximum number of iterations
    unsigned Max_iter;
 
    /// Pointer to the preconditioner
    Preconditioner* Preconditioner_pt;
 
    /// Jacobian setup time
    double Jacobian_setup_time;
 
    /// linear solver solution time
    double Solution_time;
 
    /// Preconditioner setup time
    double Preconditioner_setup_time;
 
    /// indicates whether the preconditioner should be setup before
    /// solve. Default = true;
    bool Setup_preconditioner_before_solve;
 
    /// Should we throw an error instead of just returning when we hit
    /// the max iterations?
    bool Throw_error_after_max_iter;
 
    /// Use the iterative solver as preconditioner
    bool Use_iterative_solver_as_preconditioner;
 
    /// When the iterative solver is used a preconditioner then we call
    /// the setup of solver method only once (the first time the solve
    /// method is called)
    bool First_time_solve_when_used_as_preconditioner;
  };
 
 
  ///////////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// The conjugate gradient method.
  //======================================================================
  template<typename MATRIX>
  class CG : public IterativeLinearSolver
  {
  public:
    /// Constructor
    CG()
      : Iterations(0),
        Matrix_pt(0),
        Resolving(false),
        Matrix_can_be_deleted(true)
    {
    }
 
 
    /// Destructor (cleanup storage)
    virtual ~CG()
    {
      clean_up_memory();
    }
 
    /// Broken copy constructor
    CG(const CG&) = delete;
 
    /// Broken assignment operator
    void operator=(const CG&) = delete;
 
    /// Overload disable resolve so that it cleans up memory too
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }
 
 
    /// 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 solve(Problem* const& problem_pt, DoubleVector& result);
 
    /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
    /// vector and returns the solution of the linear system.
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& solution)
    {
      // Store the matrix if required
      if ((Enable_resolve) && (!Resolving))
      {
        Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
 
        // Matrix has been passed in from the outside so we must not
        // delete it
        Matrix_can_be_deleted = false;
      }
 
      // 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(rhs.distribution_pt());
      }
 
      // Call the helper function
      this->solve_helper(matrix_pt, rhs, solution);
    }
 
    /// Re-solve the system defined by the last assembled Jacobian
    /// and the rhs vector specified here. Solution is returned in the
    /// vector result.
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    /// Number of iterations taken
    unsigned iterations() const
    {
      return Iterations;
    }
 
 
  private:
    /// General interface to solve function
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
 
    /// Cleanup data that's stored for resolve (if any has been stored)
    void clean_up_memory()
    {
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        delete Matrix_pt;
        Matrix_pt = 0;
      }
    }
 
    /// Number of iterations taken
    unsigned Iterations;
 
    /// Pointer to matrix
    MATRIX* Matrix_pt;
 
    /// Boolean flag to indicate if the solve is done in re-solve mode,
    /// bypassing setup of matrix and preconditioner
    bool Resolving;
 
    /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
    /// can be deleted.
    bool Matrix_can_be_deleted;
  };
 
 
  ///////////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// The conjugate gradient method.
  //======================================================================
  template<typename MATRIX>
  class BiCGStab : public IterativeLinearSolver
  {
  public:
    /// Constructor
    BiCGStab()
      : Iterations(0),
        Matrix_pt(0),
        Resolving(false),
        Matrix_can_be_deleted(true)
    {
    }
 
 
    /// Destructor (cleanup storage)
    virtual ~BiCGStab()
    {
      clean_up_memory();
    }
 
    /// Broken copy constructor
    BiCGStab(const BiCGStab&) = delete;
 
    /// Broken assignment operator
    void operator=(const BiCGStab&) = delete;
 
    /// Overload disable resolve so that it cleans up memory too
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }
 
    /// 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 solve(Problem* const& problem_pt, DoubleVector& result);
 
    /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
    /// vector and returns the solution of the linear system.
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& solution)
    {
      // Store the matrix if required
      if ((Enable_resolve) && (!Resolving))
      {
        Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
 
        // Matrix has been passed in from the outside so we must not
        // delete it
        Matrix_can_be_deleted = false;
      }
 
      // 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(rhs.distribution_pt());
      }
 
      // Call the helper function
      this->solve_helper(matrix_pt, rhs, solution);
    }
 
    /// Linear-algebra-type solver: Takes pointer to a matrix
    /// and rhs vector and returns the solution of the linear system
    /// Call the broken base-class version. If you want this, please
    /// implement it
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    }
 
 
    /// Re-solve the system defined by the last assembled Jacobian
    /// and the rhs vector specified here. Solution is returned in the
    /// vector result.
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
 
    /// Number of iterations taken
    unsigned iterations() const
    {
      return Iterations;
    }
 
 
  private:
    /// General interface to solve function
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    /// Cleanup data that's stored for resolve (if any has been stored)
    void clean_up_memory()
    {
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        delete Matrix_pt;
        Matrix_pt = 0;
      }
    }
 
    /// Number of iterations taken
    unsigned Iterations;
 
    /// Pointer to matrix
    MATRIX* Matrix_pt;
 
    /// Boolean flag to indicate if the solve is done in re-solve mode,
    /// bypassing setup of matrix and preconditioner
    bool Resolving;
 
    /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
    /// can be deleted.
    bool Matrix_can_be_deleted;
  };
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Smoother class:
  /// The smoother class is designed for to be used in conjunction with
  /// multigrid. The action of the smoother should reduce the high
  /// frequency errors. These methods are inefficient as stand-alone
  /// solvers.
  //====================================================================
  class Smoother : public IterativeLinearSolver
  {
  public:
    /// Empty constructor
    Smoother() : Use_as_smoother(false) {}
 
    /// Virtual empty destructor
    virtual ~Smoother() {};
 
    /// The smoother_solve function performs fixed number of iterations
    /// on the system A*result=rhs. The number of (smoothing) iterations is
    /// the same as the max. number of iterations in the underlying
    /// IterativeLinearSolver class. Note that the result vector MUST NOT
    /// re-initialised to zero (as it would typically be when the Smoother is
    /// called as a iterative linear solver).
    virtual void smoother_solve(const DoubleVector& rhs,
                                DoubleVector& result) = 0;
 
    /// Set up the smoother for the matrix specified by the pointer
    virtual void smoother_setup(DoubleMatrixBase* matrix_pt) = 0;
 
    /// Self-test to check that all the dimensions of the inputs to
    /// solve helper are consistent and everything that needs to be built, is.
    template<typename MATRIX>
    void check_validity_of_solve_helper_inputs(MATRIX* const& matrix_pt,
                                               const DoubleVector& rhs,
                                               DoubleVector& solution,
                                               const double& n_dof);
 
  protected:
    /// When a derived class object is being used as a smoother in
    /// the MG solver (or elsewhere) the residual norm does not need to be
    /// calculated because we're simply performing a fixed number of (smoothing)
    /// iterations. This boolean is used as a flag to indicate that the
    /// IterativeLinearSolver (which this class is by inheritance) is supposed
    /// to act in this way.
    bool Use_as_smoother;
  }; // End of Smoother
 
 
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
 
 
  //=========================================================================
  /// The Gauss Seidel method
  //=========================================================================
  template<typename MATRIX>
  class GS : public virtual Smoother
  {
  public:
    /// Constructor
    GS()
      : Matrix_pt(0),
        Iterations(0),
        Resolving(false),
        Matrix_can_be_deleted(true)
    {
    }
 
    /// Destructor (cleanup storage)
    virtual ~GS()
    {
      clean_up_memory();
    }
 
    /// Broken copy constructor
    GS(const GS&) = delete;
 
    /// Broken assignment operator
    void operator=(const GS&) = delete;
 
    /// Overload disable resolve so that it cleans up memory too
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    } // End of disable_resolve
 
    /// Set up the smoother for the matrix specified by the pointer
    void smoother_setup(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<MATRIX*>(matrix_pt);
    } // End of smoother_setup
 
    /// The smoother_solve function performs fixed number of iterations
    /// on the system A*result=rhs. The number of (smoothing) iterations is
    /// the same as the max. number of iterations in the underlying
    /// IterativeLinearSolver class.
    void smoother_solve(const DoubleVector& rhs, DoubleVector& result)
    {
      // If you use a smoother but you don't want to calculate the residual
      Use_as_smoother = true;
 
      // Call the helper function
      solve_helper(Matrix_pt, rhs, result);
    } // End of smoother_setup
 
    /// 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 solve(Problem* const& problem_pt, DoubleVector& result);
 
    /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
    /// vector and returns the solution of the linear system.
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& solution)
    {
      // Reset the Use_as_smoother_flag as the solver is not being used
      // as a smoother
      Use_as_smoother = false;
 
      // Set up the distribution
      this->build_distribution(rhs.distribution_pt());
 
      // Store the matrix if required
      if ((Enable_resolve) && (!Resolving))
      {
        // Upcast to the appropriate matrix type
        Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
      }
 
      // Matrix has been passed in from the outside so we must not delete it
      Matrix_can_be_deleted = false;
 
      // Call the helper function
      this->solve_helper(matrix_pt, rhs, solution);
    } // End of solve
 
    /// Linear-algebra-type solver: Takes pointer to a matrix
    /// and rhs vector and returns the solution of the linear system
    /// Call the broken base-class version. If you want this, please
    /// implement it
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    } // End of solve
 
    /// Re-solve the system defined by the last assembled Jacobian
    /// and the rhs vector specified here. Solution is returned in the
    /// vector result.
    void 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
      solve(Matrix_pt, rhs, result);
 
      // Reset re-solving flag
      Resolving = false;
    } // End of resolve
 
    /// Returns the time taken to set up the preconditioner
    double preconditioner_setup_time() const
    {
      throw OomphLibError(
        "Gauss Seidel is not a preconditionable iterative solver",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      return 0;
    } // End of preconditioner_setup_time
 
    /// Number of iterations taken
    unsigned iterations() const
    {
      return Iterations;
    } // End of iterations
 
  private:
    /// General interface to solve function
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    /// Cleanup data that's stored for resolve (if any has been stored)
    void clean_up_memory()
    {
      // If the matrix pointer isn't null and we're allowed to delete it
      // delete the matrix and assign the pointer the value NULL
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        // Destroy the matrix
        delete Matrix_pt;
 
        // Make it a null pointer
        Matrix_pt = 0;
      }
    } // End of clean_up_memory
 
    /// System matrix pointer in the format specified by the template argument
    MATRIX* Matrix_pt;
 
    /// Number of iterations taken
    unsigned Iterations;
 
    /// Boolean flag to indicate if the solve is done in re-solve mode,
    /// bypassing setup of matrix and preconditioner
    bool Resolving;
 
    /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
    /// can be deleted.
    bool Matrix_can_be_deleted;
  };
 
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
 
  //=========================================================================
  /// Explicit template specialisation of the Gauss Seidel method for
  /// compressed row format matrices
  //=========================================================================
  template<>
  class GS<CRDoubleMatrix> : public virtual Smoother
  {
  public:
    /// Constructor
    GS()
      : Matrix_pt(0),
        Iterations(0),
        Resolving(false),
        Matrix_can_be_deleted(true)
    {
    }
 
    /// Destructor (cleanup storage)
    virtual ~GS()
    {
      clean_up_memory();
    }
 
    /// Broken copy constructor
    GS(const GS&) = delete;
 
    /// Broken assignment operator
    void operator=(const GS&) = delete;
 
    /// The smoother_solve function performs fixed number of iterations
    /// on the system A*result=rhs. The number of (smoothing) iterations is
    /// the same as the max. number of iterations in the underlying
    /// IterativeLinearSolver class.
    void smoother_solve(const DoubleVector& rhs, DoubleVector& result)
    {
      // If you use a smoother but you don't want to calculate the residual
      Use_as_smoother = true;
 
      // Call the helper function
      solve_helper(Matrix_pt, rhs, result);
    } // End of smoother_solve
 
    /// Overload disable resolve so that it cleans up memory too
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    } // End of disable_resolve
 
    /// Set up the smoother for the matrix specified by the pointer
    void smoother_setup(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;
 
      // Call the generic setup helper function
      setup_helper(matrix_pt);
    } // End of smoother_setup
 
    /// Generic setup function to sort out everything that needs to be
    /// set up with regards to the input matrix
    void setup_helper(DoubleMatrixBase* matrix_pt);
 
    /// 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 solve(Problem* const& problem_pt, DoubleVector& result);
 
    /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
    /// vector and returns the solution of the linear system.
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& solution)
    {
      // Reset the Use_as_smoother_flag as the solver is not being used
      // as a smoother
      Use_as_smoother = false;
 
      // Set up the distribution
      this->build_distribution(rhs.distribution_pt());
 
      // Store the matrix if required
      if ((Enable_resolve) && (!Resolving))
      {
        // Upcast to the appropriate matrix type and sort the matrix entries
        // out so that the CRDoubleMatrix implementation of the Gauss-Seidel
        // solver can be used
        Matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
      }
      // We still need to sort the entries
      else
      {
        // 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
        dynamic_cast<CRDoubleMatrix*>(matrix_pt)->sort_entries();
 
        // Now get access to the vector Index_of_diagonal_entries
        Index_of_diagonal_entries = dynamic_cast<CRDoubleMatrix*>(matrix_pt)
                                      ->get_index_of_diagonal_entries();
      }
 
      // Matrix has been passed in from the outside so we must not delete it
      Matrix_can_be_deleted = false;
 
      // Call the helper function
      solve_helper(matrix_pt, rhs, solution);
    } // End of solve
 
    /// Linear-algebra-type solver: Takes pointer to a matrix
    /// and rhs vector and returns the solution of the linear system
    /// Call the broken base-class version. If you want this, please
    /// implement it
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    } // End of solve
 
    /// Re-solve the system defined by the last assembled Jacobian
    /// and the rhs vector specified here. Solution is returned in the
    /// vector result.
    void resolve(const DoubleVector& rhs, DoubleVector& result)
    {
      // We are re-solving
      Resolving = true;
 
#ifdef PARANOID
      // If the matrix pointer is null
      if (this->Matrix_pt == 0)
      {
        throw OomphLibError("No matrix was stored -- cannot re-solve",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Call linear algebra-style solver
      solve(Matrix_pt, rhs, result);
 
      // Reset re-solving flag
      Resolving = false;
    } // End of resolve
 
    /// Returns the time taken to set up the preconditioner
    double preconditioner_setup_time() const
    {
      // Create the error message
      std::string error_output_string = "Gauss Seidel is not a ";
      error_output_string += "preconditionable iterative solver";
 
      // Throw an error
      throw OomphLibError(
        error_output_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
 
      // Return a value so the compiler doesn't throw up an error about
      // no input being returned
      return 0;
    } // End of preconditioner_setup_time
 
    /// Number of iterations taken
    unsigned iterations() const
    {
      // Return the number of iterations
      return Iterations;
    } // End of iterations
 
  private:
    /// General interface to solve function
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    /// Clean up data that's stored for resolve (if any has been stored)
    void clean_up_memory()
    {
      // If the matrix pointer isn't null AND we're allowed to delete the
      // matrix which is only when we create the matrix ourselves
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        // Delete the matrix
        delete Matrix_pt;
 
        // Assign the associated pointer the value NULL
        Matrix_pt = 0;
      }
    } // End of clean_up_memory
 
    /// System matrix pointer in the format specified by the template argument
    CRDoubleMatrix* Matrix_pt;
 
    /// Number of iterations taken
    unsigned Iterations;
 
    /// Boolean flag to indicate if the solve is done in re-solve mode,
    /// bypassing setup of matrix and preconditioner
    bool Resolving;
 
    /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
    /// can be deleted.
    bool Matrix_can_be_deleted;
 
    /// Vector whose i'th entry contains the index of the last entry
    /// below or on the diagonal of the i'th row of the matrix
    Vector<int> Index_of_diagonal_entries;
  };
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
  //=========================================================================
  /// Damped Jacobi "solver" templated by matrix type. The "solver"
  /// exists in many different incarnations: It's an IterativeLinearSolver,
  /// and a Smoother, all of which use the same basic iteration.
  //=========================================================================
  template<typename MATRIX>
  class DampedJacobi : public virtual Smoother
  {
  public:
    /// Empty constructor
    DampedJacobi(const double& omega = 2.0 / 3.0) : Matrix_can_be_deleted(true)
    {
      // Damping factor
      Omega = omega;
    }
 
    /// Empty destructor
    ~DampedJacobi()
    {
      // Run the generic clean up function
      clean_up_memory();
    }
 
    /// Broken copy constructor
    DampedJacobi(const DampedJacobi&) = delete;
 
    /// Broken assignment operator
    void operator=(const DampedJacobi&) = delete;
 
    /// Cleanup data that's stored for resolve (if any has been stored)
    void clean_up_memory()
    {
      // If the matrix pointer isn't null AND we're allowed to delete the
      // matrix which is only when we create the matrix ourselves
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        // Delete the matrix
        delete Matrix_pt;
 
        // Assign the associated pointer the value NULL
        Matrix_pt = 0;
      }
    } // End of clean_up_memory
 
    /// Setup: Pass pointer to the matrix and store in cast form
    void smoother_setup(DoubleMatrixBase* matrix_pt)
    {
      // Assume the matrix has been passed in from the outside so we must not
      // delete it
      Matrix_can_be_deleted = false;
 
      // Upcast to the appropriate matrix type
      Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
 
      // Extract the diagonal entries of the matrix and store them
      extract_diagonal_entries(matrix_pt);
    } // End of smoother_setup
 
    /// Function to extract the diagonal entries from the matrix
    void extract_diagonal_entries(DoubleMatrixBase* matrix_pt)
    {
      // If we're using a CRDoubleMatrix object
      if (dynamic_cast<CRDoubleMatrix*>(matrix_pt))
      {
        // The matrix diagonal (we need this when we need to calculate inv(D)
        // where D is the diagonal of A and it remains the same for all uses
        // of the iterative scheme so we can store it and call it in each
        // iteration)
        Matrix_diagonal =
          dynamic_cast<CRDoubleMatrix*>(Matrix_pt)->diagonal_entries();
      }
      // If we're using a complex matrix then diagonal entries has to be a
      // complex vector rather than a vector of doubles.
      else if (dynamic_cast<CCDoubleMatrix*>(matrix_pt))
      {
        // Make an ostringstream object to create an error message
        std::ostringstream error_message_stream;
 
        // Create the error message
        error_message_stream << "Damped Jacobi can only cater to real-valued "
                             << "matrices. If you require a complex-valued "
                             << "version, please write this yourself. "
                             << "It is likely that the only difference will be "
                             << "the use of complex vectors.";
 
        // Throw an error
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      // Just extract the diagonal entries normally
      else
      {
        // Calculate the number of rows in the matrix
        unsigned n_row = Matrix_pt->nrow();
 
        // Loop over the rows of the matrix
        for (unsigned i = 0; i < n_row; i++)
        {
          // Assign the i-th value of Matrix_diagonal
          Matrix_diagonal[i] = (*Matrix_pt)(i, i);
        }
      } // if (dynamic_cast<CRDoubleMatrix*>(matrix_pt))
 
      // Calculate the n.d.o.f.
      unsigned n_dof = Matrix_diagonal.size();
 
      // Find the reciprocal of the entries of Matrix_diagonal
      for (unsigned i = 0; i < n_dof; i++)
      {
        Matrix_diagonal[i] = 1.0 / Matrix_diagonal[i];
      }
    } // End of extract_diagonal_entries
 
    /// The smoother_solve function performs fixed number of iterations
    /// on the system A*result=rhs. The number of (smoothing) iterations is
    /// the same as the max. number of iterations in the underlying
    /// IterativeLinearSolver class.
    void smoother_solve(const DoubleVector& rhs, DoubleVector& solution)
    {
      // If you use a smoother but you don't want to calculate the residual
      Use_as_smoother = true;
 
      // Call the helper function
      solve_helper(Matrix_pt, rhs, solution);
    } // End of smoother_solve
 
    /// Use damped Jacobi iteration as an IterativeLinearSolver:
    /// This obtains the Jacobian matrix J and the residual vector r
    /// (needed for the Newton method) from the problem's get_jacobian
    /// function and returns the result of Jx=r.
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
    /// vector and returns the solution of the linear system.
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& solution)
    {
      // Matrix has been passed in from the outside so we must not delete it
      Matrix_can_be_deleted = false;
 
      // Indicate that the solver is not being used as a smoother
      Use_as_smoother = false;
 
      // Set up the distribution
      this->build_distribution(rhs.distribution_pt());
 
      // Store the matrix if required
      if ((Enable_resolve) && (!Resolving))
      {
        // Upcast to the appropriate matrix type
        Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
      }
 
      // Extract the diagonal entries of the matrix and store them
      extract_diagonal_entries(matrix_pt);
 
      // Call the helper function
      solve_helper(matrix_pt, rhs, solution);
    } // End of solve
 
    /// Re-solve the system defined by the last assembled Jacobian
    /// and the rhs vector specified here. Solution is returned in the
    /// vector result.
    void 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
      solve(Matrix_pt, rhs, result);
 
      // Reset re-solving flag
      Resolving = false;
    } // End of resolve
 
    /// Number of iterations taken
    unsigned iterations() const
    {
      // Return the value of Iterations
      return Iterations;
    } // End of iterations
 
  private:
    /// This is where the actual work is done -- different
    /// implementations for different matrix types.
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    /// Pointer to the matrix
    MATRIX* Matrix_pt;
 
    /// Vector containing the diagonal entries of A
    Vector<double> Matrix_diagonal;
 
    /// Boolean flag to indicate if the solve is done in re-solve mode,
    /// bypassing setup of matrix and preconditioner
    bool Resolving;
 
    /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
    /// can be deleted.
    bool Matrix_can_be_deleted;
 
    /// Number of iterations taken
    unsigned Iterations;
 
    /// Damping factor
    double Omega;
  };
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// The GMRES method.
  //======================================================================
  template<typename MATRIX>
  class GMRES : public IterativeLinearSolver
  {
  public:
    /// Constructor
    GMRES()
      : Iterations(0),
        Matrix_pt(0),
        Resolving(false),
        Matrix_can_be_deleted(true),
        Preconditioner_application_time(0.0)
    {
      Preconditioner_LHS = true;
      Iteration_restart = false;
    }
 
    /// Destructor (cleanup storage)
    virtual ~GMRES()
    {
      clean_up_memory();
    }
 
    /// Broken copy constructor
    GMRES(const GMRES&) = delete;
 
    /// Broken assignment operator
    void operator=(const GMRES&) = delete;
 
    /// Overload disable resolve so that it cleans up memory too
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }
 
    /// function to enable the computation of the gradient
    void enable_computation_of_gradient()
    {
      Compute_gradient = true;
    }
 
    /// 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 solve(Problem* const& problem_pt, DoubleVector& result);
 
    /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
    /// vector and returns the solution of the linear system.
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& solution)
    {
      // setup the distribution
      this->build_distribution(rhs.distribution_pt());
 
      // Store the matrix if required
      if ((Enable_resolve) && (!Resolving))
      {
        Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
 
        // Matrix has been passed in from the outside so we must not
        // delete it
        Matrix_can_be_deleted = false;
      }
 
      // Call the helper function
      this->solve_helper(matrix_pt, rhs, solution);
    }
 
 
    /// Linear-algebra-type solver: Takes pointer to a matrix
    /// and rhs vector and returns the solution of the linear system
    /// Call the broken base-class version. If you want this, please
    /// implement it
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    }
 
    /// Re-solve the system defined by the last assembled Jacobian
    /// and the rhs vector specified here. Solution is returned in the
    /// vector result.
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    /// Number of iterations taken
    unsigned iterations() const
    {
      return Iterations;
    }
 
    /// access function indicating whether restarted GMRES is used
    bool iteration_restart() const
    {
      return Iteration_restart;
    }
 
    /// switches on iteration restarting and takes as an argument the
    /// number of iterations after which the construction of the
    /// orthogonalisation basis vectors should be restarted
    void enable_iteration_restart(const unsigned& restart)
    {
      Restart = restart;
      Iteration_restart = true;
    }
 
    /// switches off iteration restart
    void disable_iteration_restart()
    {
      Iteration_restart = false;
    }
 
    /// Set left preconditioning (the default)
    void set_preconditioner_LHS()
    {
      Preconditioner_LHS = true;
    }
 
    /// Enable right preconditioning
    void set_preconditioner_RHS()
    {
      Preconditioner_LHS = false;
    }
 
  private:
    /// General interface to solve function
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    /// Cleanup data that's stored for resolve (if any has been stored)
    void clean_up_memory()
    {
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        delete Matrix_pt;
        Matrix_pt = 0;
      }
    }
 
    /// Helper function to update the result vector using the result,
    /// x=x_0+V_m*y
    void update(const unsigned& k,
                const Vector<Vector<double>>& H,
                const Vector<double>& s,
                const Vector<DoubleVector>& v,
                DoubleVector& x)
    {
      // Make a local copy of s
      Vector<double> y(s);
 
      // Backsolve:
      for (int i = int(k); i >= 0; i--)
      {
        // Divide the i-th entry of y by the i-th diagonal entry of H
        y[i] /= H[i][i];
 
        // Loop over the previous values of y and update
        for (int j = i - 1; j >= 0; j--)
        {
          // Update the j-th entry of y
          y[j] -= H[i][j] * y[i];
        }
      } // for (int i=int(k);i>=0;i--)
 
      // Store the number of rows in the result vector
      unsigned n_x = x.nrow();
 
      // Build a temporary vector with entries initialised to 0.0
      DoubleVector temp(x.distribution_pt(), 0.0);
 
      // Build a temporary vector with entries initialised to 0.0
      DoubleVector z(x.distribution_pt(), 0.0);
 
      // Get access to the underlying values
      double* temp_pt = temp.values_pt();
 
      // Calculate x=Vy
      for (unsigned j = 0; j <= k; j++)
      {
        // Get access to j-th column of v
        const double* vj_pt = v[j].values_pt();
 
        // Loop over the entries of the vector, temp
        for (unsigned i = 0; i < n_x; i++)
        {
          temp_pt[i] += vj_pt[i] * y[j];
        }
      } // for (unsigned j=0;j<=k;j++)
 
      // If we're using LHS preconditioning
      if (Preconditioner_LHS)
      {
        // Since we're using LHS preconditioning the preconditioner is applied
        // to the matrix and RHS vector so we simply update the value of x
        x += temp;
      }
      // If we're using RHS preconditioning
      else
      {
        // Start the timer
        double t_start_prec = TimingHelpers::timer();
 
        // Since we're using RHS preconditioning the preconditioner is applied
        // to the solution vector
        preconditioner_pt()->preconditioner_solve(temp, z);
 
        // Calculate the time taken for the preconditioner solve
        Preconditioner_application_time +=
          (TimingHelpers::timer() - t_start_prec);
 
        // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
        // sparse linear systems", p.284]
        x += z;
      }
    } // End of update
 
    /// Helper function: Generate a plane rotation. This is done by
    /// finding the values of \f$ \cos(\theta) \f$ (i.e. cs) and
    /// \f$ \sin(\theta) \f$  (i.e. sn) such that:
    /// \f[ \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix} = \begin{bmatrix} r \newline 0 \end{bmatrix}, \f]
    /// where \f$ r=\sqrt{pow(dx,2)+pow(dy,2)} \f$. The values of a and b are
    /// given by:
    /// \f[ \cos\theta=\dfrac{dx}{\sqrt{pow(dx,2)+pow(dy,2)}}, \f]
    /// and
    /// \f[ \sin\theta=\dfrac{dy}{\sqrt{pow(dx,2)+pow(dy,2)}}. \f]
    /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.192
    void generate_plane_rotation(double& dx, double& dy, double& cs, double& sn)
    {
      // If dy=0 then we do not need to apply a rotation
      if (dy == 0.0)
      {
        // Using theta=0 gives cos(theta)=1
        cs = 1.0;
 
        // Using theta=0 gives sin(theta)=0
        sn = 0.0;
      }
      // If dx or dy is large using the normal form of calculting cs and sn
      // is naive since this may overflow or underflow so instead we calculate
      // r=sqrt(pow(dx,2)+pow(dy,2)) by using r=|dy|sqrt(1+pow(dx/dy,2)) if
      // |dy|>|dx| [see <A
      // HREF=https://en.wikipedia.org/wiki/Hypot">Hypot</A>.].
      else if (fabs(dy) > fabs(dx))
      {
        // Since |dy|>|dx| calculate the ratio dx/dy
        double temp = dx / dy;
 
        // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))
        sn = 1.0 / sqrt(1.0 + temp * temp);
 
        // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))=(dx/dy)*sin(theta)
        cs = temp * sn;
      }
      // Otherwise, we have |dx|>=|dy| so to, again, avoid overflow or underflow
      // calculate the values of cs and sn using the method above
      else
      {
        // Since |dx|>=|dy| calculate the ratio dy/dx
        double temp = dy / dx;
 
        // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))
        cs = 1.0 / sqrt(1.0 + temp * temp);
 
        // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))=(dy/dx)*cos(theta)
        sn = temp * cs;
      }
    } // End of generate_plane_rotation
 
    /// Helper function: Apply plane rotation. This is done using the
    /// update:
    /// \f[ \begin{bmatrix} dx \newline dy \end{bmatrix} \leftarrow \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix}. \f]
    void apply_plane_rotation(double& dx, double& dy, double& cs, double& sn)
    {
      // Calculate the value of dx but don't update it yet
      double temp = cs * dx + sn * dy;
 
      // Set the value of dy
      dy = -sn * dx + cs * dy;
 
      // Set the value of dx using the correct values of dx and dy
      dx = temp;
    }
 
    /// Number of iterations taken
    unsigned Iterations;
 
    /// The number of iterations before the iteration proceedure is
    /// restarted if iteration restart is used
    unsigned Restart;
 
    /// boolean indicating if iteration restarting is used
    bool Iteration_restart;
 
    /// Pointer to matrix
    MATRIX* Matrix_pt;
 
    /// Boolean flag to indicate if the solve is done in re-solve mode,
    /// bypassing setup of matrix and preconditioner
    bool Resolving;
 
    /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
    /// can be deleted.
    bool Matrix_can_be_deleted;
 
    /// boolean indicating use of left hand preconditioning (if true)
    /// or right hand preconditioning (if false)
    bool Preconditioner_LHS;
 
    /// Storage for the time spent applying the preconditioner
    double Preconditioner_application_time;
  };
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// The GMRES method.
  //======================================================================
  class AugmentedProblemGMRES : public IterativeLinearSolver
  {
  public:
    /// Constructor
    AugmentedProblemGMRES(DoubleVector* b_pt,
                          DoubleVector* c_pt,
                          double* x_pt,
                          double* rhs_pt)
      : Iterations(0),
        Matrix_pt(0),
        B_pt(b_pt),
        C_pt(c_pt),
        X_pt(x_pt),
        Rhs_pt(rhs_pt),
        Schur_complement_scalar(1.0),
        Resolving(false),
        Matrix_can_be_deleted(true)
    {
      Preconditioner_LHS = true;
      Iteration_restart = false;
    }
 
    /// Destructor: Clean up storage
    virtual ~AugmentedProblemGMRES()
    {
      clean_up_memory();
    }
 
    /// Broken copy constructor
    AugmentedProblemGMRES(const AugmentedProblemGMRES&) = delete;
 
    /// Broken assignment operator
    void operator=(const AugmentedProblemGMRES&) = delete;
 
    /// Overload disable resolve so that it cleans up memory too
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    } // End of disable_resolve
 
    /// 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 solve(Problem* const& problem_pt, DoubleVector& result);
 
    /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
    /// vector and returns the solution of the linear system.
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& solution)
    {
      // Upcast to a CRDoubleMatrix
      CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
      // If we can't upcast to a CRDoubleMatrix then this won't work...
      if (cr_matrix_pt == 0)
      {
        // Throw an error
        throw OomphLibError("Can't upcast input matrix to a CRDoubleMatrix!",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Set up the distribution
      this->build_distribution(cr_matrix_pt->distribution_pt());
 
      // Store the matrix if required
      if ((Enable_resolve) && (!Resolving))
      {
        // Store a pointer to the upcasted matrix
        Matrix_pt = cr_matrix_pt;
 
        // Matrix has been passed in from the outside so we must not
        // delete it
        Matrix_can_be_deleted = false;
      }
 
      // Call the helper function
      this->solve_helper(matrix_pt, rhs, solution);
    } // End of solve
 
    /// Linear-algebra-type solver: Takes pointer to a matrix
    /// and rhs vector and returns the solution of the linear system
    /// Call the broken base-class version. If you want this, please
    /// implement it
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    }
 
    /// Re-solve the system defined by the last assembled Jacobian
    /// and the rhs vector specified here. Solution is returned in the
    /// vector result.
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    /// Number of iterations taken
    unsigned iterations() const
    {
      return Iterations;
    }
 
    /// access function indicating whether restarted GMRES is used
    bool iteration_restart() const
    {
      return Iteration_restart;
    }
 
    /// switches on iteration restarting and takes as an argument the
    /// number of iterations after which the construction of the
    /// orthogonalisation basis vectors should be restarted
    void enable_iteration_restart(const unsigned& restart)
    {
      Restart = restart;
      Iteration_restart = true;
    }
 
    /// Switches off iteration restart
    void disable_iteration_restart()
    {
      Iteration_restart = false;
    } // End of disable_iteration_restart
 
    /// Set left preconditioning (the default)
    void set_preconditioner_LHS()
    {
      Preconditioner_LHS = true;
    }
 
    /// Enable right preconditioning
    void set_preconditioner_RHS()
    {
      Preconditioner_LHS = false;
    }
 
  private:
    /// Cleanup data that's stored for resolve (if any has been stored)
    void clean_up_memory()
    {
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        delete Matrix_pt;
        Matrix_pt = 0;
      }
    } // End of clean_up_memory
 
    /// Multiply the vector x by the augmented system matrix
    void augmented_matrix_multiply(CRDoubleMatrix* matrix_pt,
                                   const DoubleVector& x,
                                   DoubleVector& soln)
    {
      // How many dofs are there in the non-augmented system?
      unsigned n_dof = matrix_pt->nrow();
 
      // Allocate space for the non-augmented component of x
      DoubleVector x_small(this->distribution_pt(), 0.0);
 
      // Loop over the first n_dof entries of x
      for (unsigned i = 0; i < n_dof; i++)
      {
        // Copy the i-th entry over
        x_small[i] = x[i];
      }
 
      // Allocate space for the product of the matrix and x_small
      DoubleVector a_prod_xsmall(this->distribution_pt(), 0.0);
 
      // Now multiply the matrix and x_small
      matrix_pt->multiply(x_small, a_prod_xsmall);
 
      // Get the scalar component of x associated with the system augmentation
      double y = x[n_dof];
 
      // Loop over the first n_dof entries of soln
      for (unsigned i = 0; i < n_dof; i++)
      {
        // Compute the i-th entry
        soln[i] = a_prod_xsmall[i] + y * (*B_pt)[i];
      }
 
      // Compute the final entry of soln
      soln[n_dof] = C_pt->dot(x_small);
    } // End of augmented_matrix_multiply
 
    /// Apply the block-diagonal Schur complement preconditioner to
    /// compute the LHS which has size N+1 (the augmented system size)
    void apply_schur_complement_preconditioner(const DoubleVector& rhs,
                                               DoubleVector& soln)
    {
      // How many dofs are there in the non-augmented system?
      unsigned n_dof = this->distribution_pt()->nrow();
 
      // Compute the final entry of r first
      soln[n_dof] = rhs[n_dof] / Schur_complement_scalar;
 
      // Do we want to use a block-diagonal approximation? The default is to
      // use a block upper-triangular approximation for the preconditioner
      bool use_block_diagonal_preconditioner = false;
 
      // Allocate space for (the non-augmented part of) r satisfying b-Jx = Mr
      DoubleVector r_small(this->distribution_pt(), 0.0);
 
      // Allocate space for the non-augmented system part of the RHS vector
      DoubleVector rhs_small(this->distribution_pt(), 0.0);
 
      // Loop over the first n_dof entries of the RHS vector
      for (unsigned i = 0; i < n_dof; i++)
      {
        // If we want to use the block-diagonal preconditioner
        if (use_block_diagonal_preconditioner)
        {
          // Copy the i-th entry over
          rhs_small[i] = rhs[i];
        }
        // If we want to use the block upper-triangular preconditioner
        else
        {
          // Copy the i-th entry over
          rhs_small[i] = rhs[i] - soln[n_dof] * (*B_pt)[i];
        }
      } // for (unsigned i=0;i<n_dof;i++)
 
      // Apply the preconditioner
      preconditioner_pt()->preconditioner_solve(rhs_small, r_small);
 
      // Loop over the first n_dof entries of the solution vector
      for (unsigned i = 0; i < n_dof; i++)
      {
        // Copy the i-th entry over
        soln[i] = r_small[i];
      }
    } // End of apply_schur_complement_preconditioner
 
    /// General interface to solve function
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    /// Helper function to update the result vector using the result,
    /// x=x_0+V_m*y
    void update(const unsigned& k,
                const Vector<Vector<double>>& H,
                const Vector<double>& s,
                const Vector<DoubleVector>& v,
                DoubleVector& x)
    {
      // Make a local copy of s
      Vector<double> y(s);
 
      // Backsolve:
      for (int i = int(k); i >= 0; i--)
      {
        // Divide the i-th entry of y by the i-th diagonal entry of H
        y[i] /= H[i][i];
 
        // Loop over the previous values of y and update
        for (int j = i - 1; j >= 0; j--)
        {
          // Update the j-th entry of y
          y[j] -= H[i][j] * y[i];
        }
      } // for (int i=int(k);i>=0;i--)
 
      // Store the number of rows in the result vector
      unsigned n_x = x.nrow();
 
      // Build a temporary vector with entries initialised to 0.0
      DoubleVector temp(x.distribution_pt(), 0.0);
 
      // Build a temporary vector with entries initialised to 0.0
      DoubleVector z(x.distribution_pt(), 0.0);
 
      // Get access to the underlying values
      double* temp_pt = temp.values_pt();
 
      // Calculate x=Vy
      for (unsigned j = 0; j <= k; j++)
      {
        // Get access to j-th column of v
        const double* vj_pt = v[j].values_pt();
 
        // Loop over the entries of the vector, temp
        for (unsigned i = 0; i < n_x; i++)
        {
          temp_pt[i] += vj_pt[i] * y[j];
        }
      } // for (unsigned j=0;j<=k;j++)
 
      // If we're using LHS preconditioning
      if (Preconditioner_LHS)
      {
        // Since we're using LHS preconditioning the preconditioner is applied
        // to the matrix and RHS vector so we simply update the value of x
        x += temp;
      }
      // If we're using RHS preconditioning
      else
      {
        // Since we're using RHS preconditioning the preconditioner is applied
        // to the solution vector
        apply_schur_complement_preconditioner(temp, z);
 
        // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
        // sparse linear systems", p.284]
        x += z;
      }
    } // End of update
 
    /// Helper function: Generate a plane rotation. This is done by
    /// finding the values of \f$ \cos(\theta) \f$ (i.e. cs) and
    /// \f$ \sin(\theta) \f$  (i.e. sn) such that:
    /// \f[ \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix} = \begin{bmatrix} r \newline 0 \end{bmatrix}, \f]
    /// where \f$ r=\sqrt{pow(dx,2)+pow(dy,2)} \f$. The values of a and b are
    /// given by:
    /// \f[ \cos\theta=\dfrac{dx}{\sqrt{pow(dx,2)+pow(dy,2)}}, \f]
    /// and
    /// \f[ \sin\theta=\dfrac{dy}{\sqrt{pow(dx,2)+pow(dy,2)}}. \f]
    /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.192
    void generate_plane_rotation(double& dx, double& dy, double& cs, double& sn)
    {
      // If dy=0 then we do not need to apply a rotation
      if (dy == 0.0)
      {
        // Using theta=0 gives cos(theta)=1
        cs = 1.0;
 
        // Using theta=0 gives sin(theta)=0
        sn = 0.0;
      }
      // If dx or dy is large using the normal form of calculting cs and sn
      // is naive since this may overflow or underflow so instead we calculate
      // r=sqrt(pow(dx,2)+pow(dy,2)) by using r=|dy|sqrt(1+pow(dx/dy,2)) if
      // |dy|>|dx| [see <A
      // HREF=https://en.wikipedia.org/wiki/Hypot">Hypot</A>.].
      else if (fabs(dy) > fabs(dx))
      {
        // Since |dy|>|dx| calculate the ratio dx/dy
        double temp = dx / dy;
 
        // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))
        sn = 1.0 / sqrt(1.0 + temp * temp);
 
        // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))=(dx/dy)*sin(theta)
        cs = temp * sn;
      }
      // Otherwise, we have |dx|>=|dy| so to, again, avoid overflow or underflow
      // calculate the values of cs and sn using the method above
      else
      {
        // Since |dx|>=|dy| calculate the ratio dy/dx
        double temp = dy / dx;
 
        // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))
        cs = 1.0 / sqrt(1.0 + temp * temp);
 
        // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))=(dy/dx)*cos(theta)
        sn = temp * cs;
      }
    } // End of generate_plane_rotation
 
    /// Helper function: Apply plane rotation. This is done using the
    /// update:
    /// \f[ \begin{bmatrix} dx \newline dy \end{bmatrix} \leftarrow \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix}. \f]
    void apply_plane_rotation(double& dx, double& dy, double& cs, double& sn)
    {
      // Calculate the value of dx but don't update it yet
      double temp = cs * dx + sn * dy;
 
      // Set the value of dy
      dy = -sn * dx + cs * dy;
 
      // Set the value of dx using the correct values of dx and dy
      dx = temp;
    }
 
    /// Number of iterations taken
    unsigned Iterations;
 
    /// The number of iterations before the iteration proceedure is
    /// restarted if iteration restart is used
    unsigned Restart;
 
    /// boolean indicating if iteration restarting is used
    bool Iteration_restart;
 
    /// Pointer to matrix
    CRDoubleMatrix* Matrix_pt;
 
    /// Pointer to the column vector in the bordered system
    DoubleVector* B_pt;
 
    /// Pointer to the row vector in the bordered system
    DoubleVector* C_pt;
 
    /// Pointer to the last entry of the LHS vector in the bordered system
    double* X_pt;
 
    /// Pointer to the last entry of the RHS vector in the bordered system
    double* Rhs_pt;
 
    /// The scalar component of the Schur complement preconditioner
    double Schur_complement_scalar;
 
    /// Boolean flag to indicate if the solve is done in re-solve mode,
    /// bypassing setup of matrix and preconditioner
    bool Resolving;
 
    /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
    /// can be deleted.
    bool Matrix_can_be_deleted;
 
    /// boolean indicating use of left hand preconditioning (if true)
    /// or right hand preconditioning (if false)
    bool Preconditioner_LHS;
  };
} // End of namespace oomph
 
#endif