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-2022 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 generated by autoconfig
#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
{
  //===========================================================================
  //===========================================================================
  class IterativeLinearSolver : public LinearSolver
  {
  public:
    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;
    }
 
    IterativeLinearSolver(const IterativeLinearSolver&) = delete;
 
    void operator=(const IterativeLinearSolver&) = delete;
 
    virtual ~IterativeLinearSolver() {}
 
    Preconditioner*& preconditioner_pt()
    {
      return Preconditioner_pt;
    }
 
    Preconditioner* const& preconditioner_pt() const
    {
      return Preconditioner_pt;
    }
 
    double& tolerance()
    {
      return Tolerance;
    }
 
    unsigned& max_iter()
    {
      return Max_iter;
    }
 
    virtual unsigned iterations() const = 0;
 
    void enable_doc_convergence_history()
    {
      Doc_convergence_history = true;
    }
 
    void disable_doc_convergence_history()
    {
      Doc_convergence_history = false;
    }
 
    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;
    }
 
    void close_convergence_history_file_stream()
    {
      if (Output_file_stream.is_open()) Output_file_stream.close();
    }
 
    double jacobian_setup_time() const
    {
      return Jacobian_setup_time;
    }
 
    double linear_solver_solution_time() const
    {
      return Solution_time;
    }
 
    virtual double preconditioner_setup_time() const
    {
      return Preconditioner_setup_time;
    }
 
    void enable_setup_preconditioner_before_solve()
    {
      Setup_preconditioner_before_solve = true;
    }
 
    void disable_setup_preconditioner_before_solve()
    {
      Setup_preconditioner_before_solve = false;
    }
 
    void enable_error_after_max_iter()
    {
      Throw_error_after_max_iter = true;
    }
 
    void disable_error_after_max_iter()
    {
      Throw_error_after_max_iter = false;
    }
 
    void enable_iterative_solver_as_preconditioner()
    {
      Use_iterative_solver_as_preconditioner = true;
    }
 
    void disable_iterative_solver_as_preconditioner()
    {
      Use_iterative_solver_as_preconditioner = false;
    }
 
  protected:
    bool Doc_convergence_history;
 
    std::ofstream Output_file_stream;
 
    static IdentityPreconditioner Default_preconditioner;
 
    double Tolerance;
 
    unsigned Max_iter;
 
    Preconditioner* Preconditioner_pt;
 
    double Jacobian_setup_time;
 
    double Solution_time;
 
    double Preconditioner_setup_time;
 
    bool Setup_preconditioner_before_solve;
 
    bool Throw_error_after_max_iter;
 
    bool Use_iterative_solver_as_preconditioner;
 
    bool First_time_solve_when_used_as_preconditioner;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  template<typename MATRIX>
  class CG : public IterativeLinearSolver
  {
  public:
    CG()
      : Iterations(0),
        Matrix_pt(0),
        Resolving(false),
        Matrix_can_be_deleted(true)
    {
    }
 
 
    virtual ~CG()
    {
      clean_up_memory();
    }
 
    CG(const CG&) = delete;
 
    void operator=(const CG&) = delete;
 
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }
 
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    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);
    }
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    unsigned iterations() const
    {
      return Iterations;
    }
 
 
  private:
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
 
    void clean_up_memory()
    {
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        delete Matrix_pt;
        Matrix_pt = 0;
      }
    }
 
    unsigned Iterations;
 
    MATRIX* Matrix_pt;
 
    bool Resolving;
 
    bool Matrix_can_be_deleted;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  template<typename MATRIX>
  class BiCGStab : public IterativeLinearSolver
  {
  public:
    BiCGStab()
      : Iterations(0),
        Matrix_pt(0),
        Resolving(false),
        Matrix_can_be_deleted(true)
    {
    }
 
 
    virtual ~BiCGStab()
    {
      clean_up_memory();
    }
 
    BiCGStab(const BiCGStab&) = delete;
 
    void operator=(const BiCGStab&) = delete;
 
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    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);
    }
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    }
 
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
 
    unsigned iterations() const
    {
      return Iterations;
    }
 
 
  private:
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    void clean_up_memory()
    {
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        delete Matrix_pt;
        Matrix_pt = 0;
      }
    }
 
    unsigned Iterations;
 
    MATRIX* Matrix_pt;
 
    bool Resolving;
 
    bool Matrix_can_be_deleted;
  };
 
 
 
 
  //====================================================================
  //====================================================================
  class Smoother : public IterativeLinearSolver
  {
  public:
    Smoother() : Use_as_smoother(false) {}
 
    virtual ~Smoother(){};
 
    virtual void smoother_solve(const DoubleVector& rhs,
                                DoubleVector& result) = 0;
 
    virtual void smoother_setup(DoubleMatrixBase* matrix_pt) = 0;
 
    template<typename MATRIX>
    void check_validity_of_solve_helper_inputs(MATRIX* const& matrix_pt,
                                               const DoubleVector& rhs,
                                               DoubleVector& solution,
                                               const double& n_dof);
 
  protected:
    bool Use_as_smoother;
  }; // End of Smoother
 
 
 
 
  //=========================================================================
  //=========================================================================
  template<typename MATRIX>
  class GS : public virtual Smoother
  {
  public:
    GS()
      : Matrix_pt(0),
        Iterations(0),
        Resolving(false),
        Matrix_can_be_deleted(true)
    {
    }
 
    virtual ~GS()
    {
      clean_up_memory();
    }
 
    GS(const GS&) = delete;
 
    void operator=(const GS&) = delete;
 
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    } // End of disable_resolve
 
    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
 
    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
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    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
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    } // End of solve
 
    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
 
    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
 
    unsigned iterations() const
    {
      return Iterations;
    } // End of iterations
 
  private:
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    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
 
    MATRIX* Matrix_pt;
 
    unsigned Iterations;
 
    bool Resolving;
 
    bool Matrix_can_be_deleted;
  };
 
 
  //=========================================================================
  //=========================================================================
  template<>
  class GS<CRDoubleMatrix> : public virtual Smoother
  {
  public:
    GS()
      : Matrix_pt(0),
        Iterations(0),
        Resolving(false),
        Matrix_can_be_deleted(true)
    {
    }
 
    virtual ~GS()
    {
      clean_up_memory();
    }
 
    GS(const GS&) = delete;
 
    void operator=(const GS&) = delete;
 
    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
 
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    } // End of disable_resolve
 
    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
 
    void setup_helper(DoubleMatrixBase* matrix_pt);
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    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
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    } // End of solve
 
    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
 
    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
 
    unsigned iterations() const
    {
      // Return the number of iterations
      return Iterations;
    } // End of iterations
 
  private:
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    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
 
    CRDoubleMatrix* Matrix_pt;
 
    unsigned Iterations;
 
    bool Resolving;
 
    bool Matrix_can_be_deleted;
 
    Vector<int> Index_of_diagonal_entries;
  };
 
 
  //=========================================================================
  //=========================================================================
  template<typename MATRIX>
  class DampedJacobi : public virtual Smoother
  {
  public:
    DampedJacobi(const double& omega = 2.0 / 3.0) : Matrix_can_be_deleted(true)
    {
      // Damping factor
      Omega = omega;
    }
 
    ~DampedJacobi()
    {
      // Run the generic clean up function
      clean_up_memory();
    }
 
    DampedJacobi(const DampedJacobi&) = delete;
 
    void operator=(const DampedJacobi&) = delete;
 
    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
 
    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
 
    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
 
    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
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    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
 
    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
 
    unsigned iterations() const
    {
      // Return the value of Iterations
      return Iterations;
    } // End of iterations
 
  private:
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    MATRIX* Matrix_pt;
 
    Vector<double> Matrix_diagonal;
 
    bool Resolving;
 
    bool Matrix_can_be_deleted;
 
    unsigned Iterations;
 
    double Omega;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  template<typename MATRIX>
  class GMRES : public IterativeLinearSolver
  {
  public:
    GMRES()
      : Iterations(0),
        Matrix_pt(0),
        Resolving(false),
        Matrix_can_be_deleted(true),
        Preconditioner_application_time(0.0)
    {
      Preconditioner_LHS = true;
      Iteration_restart = false;
    }
 
    virtual ~GMRES()
    {
      clean_up_memory();
    }
 
    GMRES(const GMRES&) = delete;
 
    void operator=(const GMRES&) = delete;
 
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }
 
    void enable_computation_of_gradient()
    {
      Compute_gradient = true;
    }
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    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);
    }
 
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    }
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    unsigned iterations() const
    {
      return Iterations;
    }
 
    bool iteration_restart() const
    {
      return Iteration_restart;
    }
 
    void enable_iteration_restart(const unsigned& restart)
    {
      Restart = restart;
      Iteration_restart = true;
    }
 
    void disable_iteration_restart()
    {
      Iteration_restart = false;
    }
 
    void set_preconditioner_LHS()
    {
      Preconditioner_LHS = true;
    }
 
    void set_preconditioner_RHS()
    {
      Preconditioner_LHS = false;
    }
 
  private:
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    void clean_up_memory()
    {
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        delete Matrix_pt;
        Matrix_pt = 0;
      }
    }
 
    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
 
    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
 
    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;
    }
 
    unsigned Iterations;
 
    unsigned Restart;
 
    bool Iteration_restart;
 
    MATRIX* Matrix_pt;
 
    bool Resolving;
 
    bool Matrix_can_be_deleted;
 
    bool Preconditioner_LHS;
 
    double Preconditioner_application_time;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  class AugmentedProblemGMRES : public IterativeLinearSolver
  {
  public:
    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;
    }
 
    virtual ~AugmentedProblemGMRES()
    {
      clean_up_memory();
    }
 
    AugmentedProblemGMRES(const AugmentedProblemGMRES&) = delete;
 
    void operator=(const AugmentedProblemGMRES&) = delete;
 
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    } // End of disable_resolve
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    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
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    }
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    unsigned iterations() const
    {
      return Iterations;
    }
 
    bool iteration_restart() const
    {
      return Iteration_restart;
    }
 
    void enable_iteration_restart(const unsigned& restart)
    {
      Restart = restart;
      Iteration_restart = true;
    }
 
    void disable_iteration_restart()
    {
      Iteration_restart = false;
    } // End of disable_iteration_restart
 
    void set_preconditioner_LHS()
    {
      Preconditioner_LHS = true;
    }
 
    void set_preconditioner_RHS()
    {
      Preconditioner_LHS = false;
    }
 
  private:
    void clean_up_memory()
    {
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        delete Matrix_pt;
        Matrix_pt = 0;
      }
    } // End of clean_up_memory
 
    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
 
    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
 
    void solve_helper(DoubleMatrixBase* const& matrix_pt,
                      const DoubleVector& rhs,
                      DoubleVector& solution);
 
    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
 
    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
 
    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;
    }
 
    unsigned Iterations;
 
    unsigned Restart;
 
    bool Iteration_restart;
 
    CRDoubleMatrix* Matrix_pt;
 
    DoubleVector* B_pt;
 
    DoubleVector* C_pt;
 
    double* X_pt;
 
    double* Rhs_pt;
 
    double Schur_complement_scalar;
 
    bool Resolving;
 
    bool Matrix_can_be_deleted;
 
    bool Preconditioner_LHS;
  };
} // End of namespace oomph
 
#endif