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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC//
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC//
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC//
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC//
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC//
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC//
//LIC//====================================================================
//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
{

//=============================================================================
/// \short 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:

  /// \short 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&)
   {
    BrokenCopy::broken_copy("IterativeLinearSolver");
   }

  /// Broken assignment operator
  void operator=(const IterativeLinearSolver&)
   {
    BrokenCopy::broken_assign("IterativeLinearSolver");
   }

  /// 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;}

  /// \short 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();
   }

  ///  \short returns the time taken to assemble the jacobian matrix and
  /// residual vector
  double jacobian_setup_time() const
   {
    return Jacobian_setup_time;
   }

  /// \short return the time taken to solve the linear system
  double linear_solver_solution_time() const
   {
    return Solution_time;
   }

  /// \short 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:

  /// \short 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;

  /// \short 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;

  /// \short indicates whether the preconditioner should be setup before solve.
  /// Default = true;
  bool Setup_preconditioner_before_solve;

  /// \short Should we throw an error instead of just returning when we hit
  /// the max iterations?
  bool Throw_error_after_max_iter;

  /// \short 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;
 };


///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////


//======================================================================
/// \short 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&)
   {
    BrokenCopy::broken_copy("CG");
   }

  /// Broken assignment operator
  void operator=(const CG&)
   {
    BrokenCopy::broken_assign("CG");
   }

  /// Overload disable resolve so that it cleans up memory too
  void disable_resolve()
   {
    LinearSolver::disable_resolve();
    clean_up_memory();
   }


  /// \short 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);

  /// \short 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);
   }

  /// \short 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;

  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
  /// bypassing setup of matrix and preconditioner
  bool Resolving;

  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
  /// can be deleted.
  bool Matrix_can_be_deleted;
 };


///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////


//======================================================================
/// \short 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&)
   {
    BrokenCopy::broken_copy("BiCGStab");
   }

  /// Broken assignment operator
  void operator=(const BiCGStab&)
   {
    BrokenCopy::broken_assign("BiCGStab");
   }



  /// Overload disable resolve so that it cleans up memory too
  void disable_resolve()
   {
    LinearSolver::disable_resolve();
    clean_up_memory();
   }

  /// \short 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);

  /// \short 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);
   }

  /// \short 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);}



  /// \short 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;

  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
  /// bypassing setup of matrix and preconditioner
  bool Resolving;

  /// \short 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(){};

  /// \short 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;

  /// \short 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:

  /// \short 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


///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////


//=========================================================================
/// \short 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&)
   {
    BrokenCopy::broken_copy("GS");
   }

  /// Broken assignment operator
  void operator=(const GS&)
   {
    BrokenCopy::broken_assign("GS");
   }

  /// 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
  
  /// \short 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

  /// \short 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);

  /// \short 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

  /// \short 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

  /// \short 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;

  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
  /// bypassing setup of matrix and preconditioner
  bool Resolving;

  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
  /// can be deleted.
  bool Matrix_can_be_deleted;

 };

///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////

//=========================================================================
/// \short 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&)
   {
    BrokenCopy::broken_copy("GS");
   }

  /// Broken assignment operator
  void operator=(const GS&)
   {
    BrokenCopy::broken_assign("GS");
   }

  /// \short 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

  /// \short 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);

  /// \short 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);

  /// \short 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

  /// \short 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

  /// \short 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;

  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
  /// bypassing setup of matrix and preconditioner
  bool Resolving;

  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
  /// can be deleted.
  bool Matrix_can_be_deleted;

  /// \short 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&)
   {
    BrokenCopy::broken_copy("DampedJacobi");
   }

  /// Broken assignment operator
  void operator=(const DampedJacobi&)
   {
    BrokenCopy::broken_assign("DampedJacobi");
   }

  /// 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

  /// \short 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

  /// \short 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);

  /// \short 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

  /// \short 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:

  /// \short 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;

  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
  /// bypassing setup of matrix and preconditioner
  bool Resolving;

  /// \short 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;
 };


///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////


//======================================================================
/// \short 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_LHS=true;
    Iteration_restart=false;
   }

  /// Destructor (cleanup storage)
  virtual ~GMRES()
   {
    clean_up_memory();
   }

  /// Broken copy constructor
  GMRES(const GMRES&)
   {
    BrokenCopy::broken_copy("GMRES");
   }

  /// Broken assignment operator
  void operator=(const GMRES&)
   {
    BrokenCopy::broken_assign("GMRES");
   }

  /// Overload disable resolve so that it cleans up memory too
  void disable_resolve()
   {
    LinearSolver::disable_resolve();
    clean_up_memory();
   }

  /// \short 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);

  /// \short 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);
   }


  /// \short 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);}

  /// \short 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;
   }

  /// \short access function indicating whether restarted GMRES is used
  bool iteration_restart() const
   {
    return Iteration_restart;
   }

  /// \short 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;
   }

  /// \short Set left preconditioning (the default)
  void set_preconditioner_LHS() {Preconditioner_LHS=true;}

  /// \short 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
    {
     // Since we're using RHS preconditioning the preconditioner is applied
     // to the solution vector
     preconditioner_pt()->preconditioner_solve(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

  /// \short Helper function: Generate a plane rotation. This is done by
  /// finding the values of \f$ \cos(\theta) \f$ (i.e. cs) and \sin(\theta)
  /// (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

  /// \short 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;

  /// \short 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;

  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
  /// bypassing setup of matrix and preconditioner
  bool Resolving;

  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
  /// can be deleted.
  bool Matrix_can_be_deleted;

  /// \short boolean indicating use of left hand preconditioning (if true)
  /// or right hand preconditioning (if false)
  bool Preconditioner_LHS;
 };
} // End of namespace oomph

#endif