iterative_linear_solver.cc

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented,
//LIC// multi-physics finite-element library, available
//LIC// at http://www.oomph-lib.org.
//LIC//
//LIC//    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//====================================================================
//The actual solve functions for Iterative solvers.


// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif

// Oomph-lib includes
#include "iterative_linear_solver.h"

// Required to force_ get templated builds of iterative solvers for
// sumofmatrices class.
#include "sum_of_matrices.h"


namespace oomph
{


//==================================================================
/// \short Default preconditioner for iterative solvers: The base
/// class for preconditioners is a fully functional (if trivial!)
/// preconditioner.
//==================================================================
 IdentityPreconditioner IterativeLinearSolver::Default_preconditioner;


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


//==================================================================
/// \short Re-solve the system defined by the last assembled Jacobian
/// and the rhs vector specified here. Solution is returned in
/// the vector result.
//==================================================================
 template<typename MATRIX>
 void BiCGStab<MATRIX>::resolve(const DoubleVector &rhs,
                                DoubleVector &result)
 {

  // We are re-solving
  Resolving=true;

#ifdef PARANOID
  if (Matrix_pt==0)
  {
   throw OomphLibError(
    "No matrix was stored -- cannot re-solve",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
#endif

  // Call linear algebra-style solver
  this->solve(Matrix_pt,rhs,result);

  // Reset re-solving flag
  Resolving=false;

 }

//==================================================================
/// Solver: Takes pointer to problem and returns the results vector
/// which contains the solution of the linear system defined by
/// the problem's fully assembled Jacobian and residual vector.
//==================================================================
 template<typename MATRIX>
 void BiCGStab<MATRIX>::solve(Problem* const &problem_pt,
                              DoubleVector &result)
 {
  // Initialise timer
#ifdef OOMPH_HAS_MPI
  double t_start = TimingHelpers::timer();
#else
  clock_t t_start = clock();
#endif

  // We're not re-solving
  Resolving=false;

  // Get rid of any previously stored data
  clean_up_memory();

  // Get Jacobian matrix in format specified by template parameter
  // and nonlinear residual vector
  Matrix_pt=new MATRIX;
  DoubleVector f;
  problem_pt->get_jacobian(f,*Matrix_pt);


  // We've made the matrix, we can delete it...
  Matrix_can_be_deleted=true;

  // Doc time for setup
#ifdef OOMPH_HAS_MPI
  double t_end = TimingHelpers::timer();
  Jacobian_setup_time= t_end-t_start;
#else
  clock_t t_end = clock();
  Jacobian_setup_time=double(t_end-t_start)/CLOCKS_PER_SEC;
#endif

  if(Doc_time)
  {
   oomph_info << "Time for setup of Jacobian [sec]: "
              << Jacobian_setup_time << std::endl;
  }

  // set the distribution
  if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt))
  {
   // the solver has the same distribution as the matrix if possible
   this->build_distribution(dynamic_cast<DistributableLinearAlgebraObject*>
                            (Matrix_pt)->distribution_pt());
  }
  else
  {
   // the solver has the same distribution as the RHS
   this->build_distribution(f.distribution_pt());
  }

  // Call linear algebra-style solver
  if((result.built()) &&
     (!(*result.distribution_pt() == *this->distribution_pt())))
  {
   LinearAlgebraDistribution
    temp_global_dist(result.distribution_pt());
   result.build(this->distribution_pt(),0.0);
   this->solve_helper(Matrix_pt,f,result);
   result.redistribute(&temp_global_dist);
  }
  else
  {
   this->solve_helper(Matrix_pt,f,result);
  }

  // Kill matrix unless it's still required for resolve
  if (!Enable_resolve) clean_up_memory();

 };




//==================================================================
/// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
/// and returns the solution of the linear system.
/// Algorithm and variable names based on "Numerical Linear Algebra
/// for High-Performance Computers" by Dongarra, Duff, Sorensen & van
/// der Vorst. SIAM  (1998), page 185.
//==================================================================
 template<typename MATRIX>
 void BiCGStab<MATRIX>::solve_helper(DoubleMatrixBase* const &matrix_pt,
                                     const DoubleVector &rhs,
                                     DoubleVector &solution)
 {
#ifdef PARANOID
  // check that the rhs vector is setup
  if (!rhs.built())
  {
   std::ostringstream error_message_stream;
   error_message_stream
    << "The vectors rhs must be setup";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

  // check that the matrix is square
  if (matrix_pt->nrow() != matrix_pt->ncol())
  {
   std::ostringstream error_message_stream;
   error_message_stream
    << "The matrix at matrix_pt must be square.";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

  // check that the matrix and the rhs vector have the same nrow()
  if (matrix_pt->nrow() != rhs.nrow())
  {
   std::ostringstream error_message_stream;
   error_message_stream
    << "The matrix and the rhs vector must have the same number of rows.";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

  // if the matrix is distributable then it too should have the same
  // communicator as the rhs vector
  DistributableLinearAlgebraObject* dist_matrix_pt =
   dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
  if (dist_matrix_pt != 0)
  {
   if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
   {
    std::ostringstream error_message_stream;
    error_message_stream
     << "The matrix matrix_pt must have the same communicator as the vectors"
     << " rhs and result must have the same communicator";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  }
  // if the matrix is not distributable then it the rhs vector should not be
  // distributed
  else
  {
   if (rhs.distribution_pt()->distributed())
   {
    std::ostringstream error_message_stream;
    error_message_stream
     << "The matrix (matrix_pt) is not distributable and therefore the rhs"
     << " vector must not be distributed";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  }
  // if the result vector is setup then check it has the same distribution
  // as the rhs
  if (solution.built())
  {
   if (!(*solution.distribution_pt() == *rhs.distribution_pt()))
   {
    std::ostringstream error_message_stream;
    error_message_stream
     << "The solution vector distribution has been setup; it must have the "
     << "same distribution as the rhs vector.";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  }
#endif

  // setup the solution if it is not
  if (!solution.distribution_pt()->built())
  {
   solution.build(this->distribution_pt(),0.0);
  }
  // zero
  else
  {
   solution.initialise(0.0);
  }

  // Get number of dofs
// unsigned n_dof=rhs.size();
  unsigned nrow_local = this->nrow_local();

  // Time solver
#ifdef OOMPH_HAS_MPI
  double t_start = TimingHelpers::timer();
#else
  clock_t t_start = clock();
#endif

  // Initialise: Zero initial guess so the initial residual is
  // equal to the RHS, i.e. the nonlinear residual
  DoubleVector residual(rhs);
  double residual_norm = residual.norm();
  double rhs_norm = residual_norm;
  if (rhs_norm==0.0) rhs_norm=1.0;
  DoubleVector x(rhs.distribution_pt(),0.0);

  // Hat residual by copy operation
  DoubleVector r_hat(residual);

  // Normalised residual
  double normalised_residual_norm=residual_norm/rhs_norm;

  // if required will document convergence history to screen or file (if
  // stream open)
  if (Doc_convergence_history)
  {
   if (!Output_file_stream.is_open())
   {
    oomph_info << 0 << " "
               << normalised_residual_norm << std::endl;
   }
   else
   {
    Output_file_stream << 0 << " "
                       << normalised_residual_norm << std::endl;
   }
  }

  // Check immediate convergence
  if (normalised_residual_norm<Tolerance)
  {
   if(Doc_time)
   {
    oomph_info << "BiCGStab converged immediately" << std::endl;
   }
   solution=x;

   // Doc time for solver
   double t_end = TimingHelpers::timer();
   Solution_time = t_end-t_start;

   if(Doc_time)
   {
    oomph_info << "Time for solve with BiCGStab  [sec]: "
               << Solution_time << std::endl;
   }
   return;
  }

  // Setup preconditioner only if we're not re-solving
  if (!Resolving)
  {
   // only setup the preconditioner if required
   if (Setup_preconditioner_before_solve)
   {
    //Setup preconditioner from the Jacobian matrix
    double t_start_prec = TimingHelpers::timer();

    preconditioner_pt()->setup(matrix_pt);

    // Doc time for setup of preconditioner
    double t_end_prec = TimingHelpers::timer();
    Preconditioner_setup_time = t_end_prec-t_start_prec;

    if(Doc_time)
    {
     oomph_info << "Time for setup of preconditioner  [sec]: "
                << Preconditioner_setup_time << std::endl;
    }
   }
  }
  else
  {
   if(Doc_time)
   {
    oomph_info << "Setup of preconditioner is bypassed in resolve mode"
               << std::endl;
   }
  }

  // Some auxiliary variables
  double rho_prev=0.0;
  double alpha=0.0;
  double omega=0.0;
  double rho,beta,dot_prod,dot_prod_tt,dot_prod_ts;
  double s_norm,r_norm;

  // Some vectors
  DoubleVector p(this->distribution_pt(),0.0),p_hat(this->distribution_pt(),0.0),
   v(this->distribution_pt(),0.0),z(this->distribution_pt(),0.0),t(this->distribution_pt(),0.0),
   s(this->distribution_pt(),0.0);

  // Loop over max. number of iterations
  for (unsigned iter=1;iter<=Max_iter;iter++)
  {
   // Dot product for rho
   rho = r_hat.dot(residual);

   // Breakdown?
   if (rho==0.0)
   {
    oomph_info << "BiCGStab has broken down after " << iter
               << " iterations" << std::endl;
    oomph_info << "Returning with current normalised residual of "
               << normalised_residual_norm << std::endl;
   }

   // First step is different
   if (iter==1)
   {
    p=residual;
   }
   else
   {
    beta=(rho/rho_prev)*(alpha/omega);
    for (unsigned i=0;i<nrow_local;i++)
    {
     p[i] = residual[i] + beta*(p[i]-omega*v[i]);
    }
   }

   // Apply precondtitioner: p_hat=P^-1*p
   preconditioner_pt()->preconditioner_solve(p,p_hat);

   // Matrix vector product: v=A*p_hat
   matrix_pt->multiply(p_hat,v);
   dot_prod = r_hat.dot(v);
   alpha=rho/dot_prod;
   for (unsigned i=0;i<nrow_local;i++)
   {
    s[i]=residual[i]-alpha*v[i];
   }
   s_norm = s.norm();

   // Normalised residual
   normalised_residual_norm=s_norm/rhs_norm;

   // if required will document convergence history to screen or file (if
   // stream open)
   if (Doc_convergence_history)
   {
    if (!Output_file_stream.is_open())
    {
     oomph_info << double(iter-0.5) << " "
                << normalised_residual_norm << std::endl;
    }
    else
    {
     Output_file_stream << double(iter-0.5) << " "
                        << normalised_residual_norm <<std::endl;
    }
   }

   // Converged?
   if (normalised_residual_norm<Tolerance)
   {
    for (unsigned i=0;i<nrow_local;i++)
    {
     solution[i]=x[i]+alpha*p_hat[i];
    }

    if(Doc_time)
    {
     oomph_info << std::endl;
     oomph_info << "BiCGStab converged. Normalised residual norm: "
                << normalised_residual_norm << std::endl;
     oomph_info << "Number of iterations to convergence: "
                << iter << std::endl;
     oomph_info << std::endl;
    }

    // Store number of iterations taken
    Iterations = iter;

    // Doc time for solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end-t_start;

    if(Doc_time)
    {
     oomph_info << "Time for solve with BiCGStab  [sec]: "
                << Solution_time << std::endl;
    }

    return;
   }

   //Apply precondtitioner: z=P^-1*s
   preconditioner_pt()->preconditioner_solve(s,z);

   // Matrix vector product: t=A*z
   matrix_pt->multiply(z,t);
   dot_prod_ts=t.dot(s);
   dot_prod_tt=t.dot(t);
   omega=dot_prod_ts/dot_prod_tt;
   for (unsigned i=0;i<nrow_local;i++)
   {
    x[i]+=alpha*p_hat[i]+omega*z[i];
    residual[i]=s[i]-omega*t[i];
   }
   r_norm = residual.norm();
   rho_prev=rho;

   // Check convergence again
   normalised_residual_norm=r_norm/rhs_norm;

   // if required will document convergence history to screen or file (if
   // stream open)
   if (Doc_convergence_history)
   {
    if (!Output_file_stream.is_open())
    {
     oomph_info <<  iter << " "
                << normalised_residual_norm << std::endl;
    }
    else
    {
     Output_file_stream << iter << " "
                        << normalised_residual_norm << std::endl;
    }
   }


   if (normalised_residual_norm<Tolerance)
   {
    if(Doc_time)
    {
     oomph_info << std::endl;
     oomph_info << "BiCGStab converged. Normalised residual norm: "
                << normalised_residual_norm << std::endl;
     oomph_info << "Number of iterations to convergence: "
                << iter << std::endl;
     oomph_info << std::endl;
    }
    solution=x;

    // Store the number of itertions taken.
    Iterations = iter;

    // Doc time for solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end-t_start;

    if(Doc_time)
    {
     oomph_info << "Time for solve with BiCGStab  [sec]: "
                << Solution_time << std::endl;
    }
    return;
   }


   // Breakdown: Omega has to be >0 for to be able to continue
   if (omega == 0.0)
   {
    oomph_info << std::endl;
    oomph_info << "BiCGStab breakdown with omega=0.0. "
               << "Normalised residual norm: "
               << normalised_residual_norm << std::endl;
    oomph_info << "Number of iterations so far: " << iter << std::endl;
    oomph_info << std::endl;
    solution=x;

    // Store the number of itertions taken.
    Iterations = iter;

    // Doc time for solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end-t_start;

    if(Doc_time)
    {
     oomph_info << "Time for solve with BiCGStab  [sec]: "
                << Solution_time << std::endl;
    }
    return;
   }



  } // end of iteration loop


  // No convergence
  oomph_info << std::endl;
  oomph_info << "BiCGStab did not converge to required tolerance! "
             << std::endl;
  oomph_info << "Returning with normalised residual norm: "
             << normalised_residual_norm << std::endl;
  oomph_info << "after " << Max_iter<< " iterations." << std::endl;
  oomph_info << std::endl;

  solution=x;

  // Doc time for solver
  double t_end = TimingHelpers::timer();
  Solution_time = t_end-t_start;

  if(Doc_time)
  {
   oomph_info << "Time for solve with BiCGStab  [sec]: "
              << Solution_time << std::endl;
  }

  if(Throw_error_after_max_iter)
  {
   std::string err = "Solver failed to converge and you requested an error";
   err += " on convergence failures.";
   throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                       OOMPH_CURRENT_FUNCTION);
  }

 }


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



//==================================================================
/// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
/// and returns the solution of the linear system.
/// Algorithm and variable names based on "Matrix Computations,
/// 2nd Ed." Golub & van Loan, John Hopkins University Press(1989),
/// page 529.
//==================================================================
 template<typename MATRIX>
 void CG<MATRIX>::solve_helper(DoubleMatrixBase* const &matrix_pt,
                               const DoubleVector& rhs,
                               DoubleVector& solution)
 {
#ifdef PARANOID
  // check that the rhs vector is setup
  if (!rhs.built())
  {
   std::ostringstream error_message_stream;
   error_message_stream
    << "The vectors rhs must be setup";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

  // check that the matrix is square
  if (matrix_pt->nrow() != matrix_pt->ncol())
  {
   std::ostringstream error_message_stream;
   error_message_stream
    << "The matrix at matrix_pt must be square.";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

  // check that the matrix and the rhs vector have the same nrow()
  if (matrix_pt->nrow() != rhs.nrow())
  {
   std::ostringstream error_message_stream;
   error_message_stream
    << "The matrix and the rhs vector must have the same number of rows.";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

  // if the matrix is distributable then it too should have the same
  // communicator as the rhs vector
  DistributableLinearAlgebraObject* dist_matrix_pt =
   dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
  if (dist_matrix_pt != 0)
  {
   if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
   {
    std::ostringstream error_message_stream;
    error_message_stream
     << "The matrix matrix_pt must have the same communicator as the vectors"
     << " rhs and result must have the same communicator";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  }
  // if the matrix is not distributable then it the rhs vector should not be
  // distributed
  else
  {
   if (rhs.distribution_pt()->distributed())
   {
    std::ostringstream error_message_stream;
    error_message_stream
     << "The matrix (matrix_pt) is not distributable and therefore the rhs"
     << " vector must not be distributed";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  }
  // if the result vector is setup then check it has the same distribution
  // as the rhs
  if (solution.built())
  {
   if (!(*solution.distribution_pt() == *rhs.distribution_pt()))
   {
    std::ostringstream error_message_stream;
    error_message_stream
     << "The solution vector distribution has been setup; it must have the "
     << "same distribution as the rhs vector.";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  }
#endif

  // setup the solution if it is not
  if (!solution.distribution_pt()->built())
  {
   solution.build(this->distribution_pt(),0.0);
  }
  // zero
  else
  {
   solution.initialise(0.0);
  }

  // Get number of dofs
// unsigned n_dof=rhs.size();
  unsigned nrow_local = this->nrow_local();

  // Initialise counter
  unsigned counter = 0;

  // Time solver
  double t_start = TimingHelpers::timer();

  // Initialise: Zero initial guess so the initial residual is
  // equal to the RHS
  DoubleVector x(this->distribution_pt(),0.0);
  DoubleVector residual(rhs);
  double residual_norm = residual.norm();
  double rhs_norm=residual_norm;
  if (rhs_norm==0.0) rhs_norm=1.0;

  // Normalised residual
  double normalised_residual_norm=residual_norm/rhs_norm;

  // if required will document convergence history to screen or file (if
  // stream open)
  if (Doc_convergence_history)
  {
   if (!Output_file_stream.is_open())
   {
    oomph_info  << 0 << " "
                << normalised_residual_norm <<std::endl;
   }
   else
   {
    Output_file_stream << 0 << " "
                       << normalised_residual_norm <<std::endl;
   }
  }

  // Check immediate convergence
  if (normalised_residual_norm<Tolerance)
  {
   if(Doc_time)
   {
    oomph_info << "CG converged immediately" << std::endl;
   }
   solution=x;

   // Doc time for solver
   double t_end = TimingHelpers::timer();
   Solution_time = t_end-t_start;

   if(Doc_time)
   {
    oomph_info << "Time for solve with CG  [sec]: "
               << Solution_time << std::endl;
   }
   return;
  }


  // Setup preconditioner only if we're not re-solving
  if (!Resolving)
  {
   // only setup the preconditioner if required
   if (Setup_preconditioner_before_solve)
   {
    //Setup preconditioner from the Jacobian matrix
    double t_start_prec = TimingHelpers::timer();

    preconditioner_pt()->setup(matrix_pt);

    // Doc time for setup of preconditioner
    double t_end_prec = TimingHelpers::timer();
    Preconditioner_setup_time = t_end_prec-t_start_prec;

    if(Doc_time)
    {
     oomph_info << "Time for setup of preconditioner  [sec]: "
                << Preconditioner_setup_time << std::endl;
    }
   }
  }
  else
  {
   if(Doc_time)
   {
    oomph_info << "Setup of preconditioner is bypassed in resolve mode"
               << std::endl;
   }
  }


  // Auxiliary vectors
// Vector<double> z(n_dof),p(n_dof),jacobian_times_p(n_dof,0.0);
  DoubleVector z(this->distribution_pt(),0.0), p(this->distribution_pt(),0.0),
   jacobian_times_p(this->distribution_pt(),0.0);

  // Auxiliary values
  double alpha,beta,rz;
  double prev_rz=0.0;

  // Main iteration
  while((normalised_residual_norm>Tolerance)&&(counter!=Max_iter))
  {

   //Apply precondtitioner: z=P^-1*r
   preconditioner_pt()->preconditioner_solve(residual,z);

   // P vector is computed differently for first and subsequent steps
   if(counter==0)
   {
    p=z;
    rz=residual.dot(z);
   }
   // Subsequent steps
   else
   {
    rz=residual.dot(z);
    beta=rz/prev_rz;
    for (unsigned i=0;i<nrow_local;i++)
    {
     p[i]=z[i]+beta*p[i];
    }
   }


   // Matrix vector product
   matrix_pt->multiply(p,jacobian_times_p);
   double pq = p.dot(jacobian_times_p);
   alpha=rz/pq;

   // Update
   prev_rz=rz;
   for(unsigned i=0;i<nrow_local;i++)
   {
    x[i]+=alpha*p[i];
    residual[i]-=alpha*jacobian_times_p[i];
   }


   //Calculate the 2norm
   residual_norm = residual.norm();

   //Difference between the initial and current 2norm residual
   normalised_residual_norm=residual_norm/rhs_norm;


   // if required will document convergence history to screen or file (if
   // stream open)
   if (Doc_convergence_history)
   {
    if (!Output_file_stream.is_open())
    {
     oomph_info <<  counter << " "
                << normalised_residual_norm << std::endl;
    }
    else
    {
     Output_file_stream << counter << " "
                        << normalised_residual_norm << std::endl;
    }
   }

   counter=counter+1;

  }//end while


  if (counter >= Max_iter)
  {
   oomph_info << std::endl;
   oomph_info << "CG did not converge to required tolerance! " << std::endl;
   oomph_info << "Returning with normalised residual norm: "
              << normalised_residual_norm << std::endl;
   oomph_info << "after " << counter << " iterations." << std::endl;
   oomph_info << std::endl;
  }
  else
  {
   if(Doc_time)
   {
    oomph_info << std::endl;
    oomph_info << "CG converged. Normalised residual norm: "
               << normalised_residual_norm << std::endl;
    oomph_info << "Number of iterations to convergence: "
               << counter << std::endl;
    oomph_info << std::endl;
   }
  }


  // Store number if iterations taken
  Iterations = counter;

  // Copy result back
  solution=x;

  // Doc time for solver
  double t_end = TimingHelpers::timer();
  Solution_time = t_end-t_start;

  if(Doc_time)
  {
   oomph_info << "Time for solve with CG  [sec]: "
              << Solution_time << std::endl;
  }

  if((counter >= Max_iter) && (Throw_error_after_max_iter))
  {
   std::string err = "Solver failed to converge and you requested an error";
   err += " on convergence failures.";
   throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                       OOMPH_CURRENT_FUNCTION);
  }

 }//end CG


//==================================================================
/// \short Re-solve the system defined by the last assembled Jacobian
/// and the rhs vector specified here. Solution is returned in
/// the vector result.
//==================================================================
 template<typename MATRIX>
 void CG<MATRIX>::resolve(const DoubleVector &rhs,
                          DoubleVector &result)
 {

  // We are re-solving
  Resolving=true;

#ifdef PARANOID
  if (Matrix_pt==0)
  {
   throw OomphLibError(
    "No matrix was stored -- cannot re-solve",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
#endif

  // Call linear algebra-style solver
  this->solve(Matrix_pt,rhs,result);

  // Reset re-solving flag
  Resolving=false;

 }

 
//==================================================================
/// Solver: Takes pointer to problem and returns the results vector
/// which contains the solution of the linear system defined by
/// the problem's fully assembled Jacobian and residual vector.
//==================================================================
 template<typename MATRIX>
 void CG<MATRIX>::solve(Problem* const &problem_pt, DoubleVector &result)
 {
  // Initialise timer
  double t_start = TimingHelpers::timer();

  // We're not re-solving
  Resolving=false;

  // Get rid of any previously stored data
  clean_up_memory();

  // Get Jacobian matrix in format specified by template parameter
  // and nonlinear residual vector
  Matrix_pt=new MATRIX;
  DoubleVector f;
  problem_pt->get_jacobian(f,*Matrix_pt);

  // We've made the matrix, we can delete it...
  Matrix_can_be_deleted=true;

  // Doc time for setup
  double t_end = TimingHelpers::timer();
  Jacobian_setup_time= t_end-t_start;

  if(Doc_time)
  {
   oomph_info << "Time for setup of Jacobian [sec]: "
              << Jacobian_setup_time << std::endl;
  }

  // set the distribution
  if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt))
  {
   // the solver has the same distribution as the matrix if possible
   this->build_distribution(dynamic_cast<DistributableLinearAlgebraObject*>
                            (Matrix_pt)->distribution_pt());
  }
  else
  {
   // the solver has the same distribution as the RHS
   this->build_distribution(f.distribution_pt());
  }

  // if the result vector is not setup
  if (!result.distribution_pt()->built())
  {
   result.build(this->distribution_pt(),0.0);
  }

  // Call linear algebra-style solver
  if (!(*result.distribution_pt() == *this->distribution_pt()))
  {
   LinearAlgebraDistribution
    temp_global_dist(result.distribution_pt());
   result.build(this->distribution_pt(),0.0);
   this->solve_helper(Matrix_pt,f,result);
   result.redistribute(&temp_global_dist);
  }
  else
  {
   this->solve_helper(Matrix_pt,f,result);
  }

  // Kill matrix unless it's still required for resolve
  if (!Enable_resolve) clean_up_memory();

 };


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


//==================================================================
/// \short Self-test to be called inside solve_helper to ensure
/// that all inputs are consistent and everything that needs to
/// be built, is.
//==================================================================
 template<typename MATRIX>
 void Smoother::check_validity_of_solve_helper_inputs(MATRIX* const &matrix_pt,
                                                      const DoubleVector& rhs,
                                                      DoubleVector& solution,
                                                      const double& n_dof)
 {
  // Check that if the matrix is not distributable then it should not
  // be distributed
  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
  {
   if (dynamic_cast<DistributableLinearAlgebraObject*>
       (matrix_pt)->distributed())
   {
    std::ostringstream error_message_stream;
    error_message_stream << "The matrix must not be distributed.";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  } // if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)!=0)
  // Check that this rhs distribution is setup
  if (!rhs.built())
  {
   std::ostringstream error_message_stream;
   error_message_stream << "The rhs vector distribution must be setup.";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
  // Check that the rhs has the right number of global rows
  if(rhs.nrow()!=n_dof)
  {
   std::ostringstream error_message_stream;
   error_message_stream << "RHS does not have the same dimension as "
                        << "the linear system";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
  // Check that the rhs is not distributed
  if (rhs.distribution_pt()->distributed())
  {
   std::ostringstream error_message_stream;
   error_message_stream << "The rhs vector must not be distributed.";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
  // Check that if the result is setup it matches the distribution
  // of the rhs
  if (solution.built())
  {
   if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
   {
    std::ostringstream error_message_stream;
    error_message_stream << "If the result distribution is setup then it "
                         << "must be the same as the rhs distribution";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  } // if (solution.built())
 } // End of check_validity_of_solve_helper_inputs

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

 
//==================================================================
/// \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.
//==================================================================
 template<typename MATRIX>
 void GS<MATRIX>::solve(Problem* const &problem_pt,DoubleVector &result)
 {
  // Reset the Use_as_smoother_flag as the solver is not being used
  // as a smoother
  Use_as_smoother=false;
  
  // Find the # of degrees of freedom (variables)
  unsigned n_dof=problem_pt->ndof();

  // Initialise timer
  double t_start=TimingHelpers::timer();

  // We're not re-solving
  Resolving=false;

  // Get rid of any previously stored data
  clean_up_memory();

  // Set up the distribution
  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),
                                 n_dof,false);

  // Assign the distribution to the LinearSolver
  this->build_distribution(dist);

  // Allocate space for the Jacobian matrix in format specified
  // by template parameter
  Matrix_pt=new MATRIX;

  // Get the nonlinear residuals vector
  DoubleVector f;

  // Assign the Jacobian and the residuals vector
  problem_pt->get_jacobian(f,*Matrix_pt);

  // We've made the matrix, we can delete it...
  Matrix_can_be_deleted=true;

  // Doc time for setup
  double t_end=TimingHelpers::timer();
  Jacobian_setup_time=t_end-t_start;

  // If time documentation is enabled
  if(Doc_time)
  {
   oomph_info << "Time for setup of Jacobian [sec]: "
              << Jacobian_setup_time << std::endl;
  }

  // Call linear algebra-style solver
  this->solve_helper(Matrix_pt,f,result);

  // Kill matrix unless it's still required for resolve
  if (!Enable_resolve) clean_up_memory();
 } // End of solve

//==================================================================
/// \short Linear-algebra-type solver: Takes pointer to a matrix and
/// rhs vector and returns the solution of the linear system.
//==================================================================
 template<typename MATRIX>
 void GS<MATRIX>::solve_helper(DoubleMatrixBase* const &matrix_pt,
                               const DoubleVector& rhs,
                               DoubleVector& solution)
 {
  // Get number of dofs
  unsigned n_dof=rhs.nrow();

#ifdef PARANOID
  // Upcast the matrix to the appropriate type
  MATRIX* tmp_matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
  
  // PARANOID Run the self-tests to check the inputs are correct
  this->check_validity_of_solve_helper_inputs<MATRIX>(tmp_matrix_pt,rhs,
                                                      solution,n_dof);

  // We don't need the pointer any more but we do still need the matrix
  // so just make tmp_matrix_pt a null pointer
  tmp_matrix_pt=0;
#endif

  // Set up the solution if it is not
  if (!solution.distribution_pt()->built())
  {
   // Build it!
   solution.build(this->distribution_pt(),0.0);
  }
  // If the solution is already set up
  else
  {
   // If we're inside the multigrid solver then as we traverse up the
   // hierarchy we use the smoother on the updated approximate solution.
   // As such, we should ONLY be resetting all the values to zero if
   // we're NOT inside the multigrid solver
   if (!Use_as_smoother)
   {
    // Initialise the vector with all entries set to zero
    solution.initialise(0.0);
   }
  } // if (!solution.distribution_pt()->built())

  // Initialise timer
  double t_start=TimingHelpers::timer();

  // Copy the solution vector into x
  DoubleVector x(solution);

  // Create a vector to hold the residual. This will only be built if
  // we're not inside the multigrid solver
  DoubleVector current_residual;

  // Variable to hold the current residual norm. Only used if we're
  // not inside the multigrid solver
  double norm_res=0.0;

  // Variables to hold the initial residual norm. Only used if we're
  // not inside the multigrid solver
  double norm_f=0.0;

  // Initialise the value of Iterations
  Iterations=0;

  // Calculate the residual only if we're not inside the multigrid solver
  if (!Use_as_smoother)
  {
   // Build the residual vector
   current_residual.build(this->distribution_pt(),0.0);

   // Calculate the residual r=b-Ax
   matrix_pt->residual(x,rhs,current_residual);

   // Calculate the 2-norm of the residual vector
   norm_res=current_residual.norm();

   // Store the initial norm
   norm_f=norm_res;
   
   // If required, document the convergence history to screen or file (if
   // the output stream open)
   if (Doc_convergence_history)
   {
    // If the output file stream isn't open
    if (!Output_file_stream.is_open())
    {
     // Output the result to screen
     oomph_info << Iterations << " " << norm_res << std::endl;
    }
    // If the output file stream is open
    else
    {
     // Document the result to file
     Output_file_stream << Iterations << " " << norm_res <<std::endl;
    }
   } // if (Doc_convergence_history)
  } // if (!Use_as_smoother)

  // Start of the main GS loop: run up to Max_iter times
  for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
  {
   // Loop over rows
   for(unsigned i=0;i<n_dof;i++)
   {
    double dummy=rhs[i];
    for (unsigned j=0;j<i;j++)
    {
     dummy-=(*(matrix_pt))(i,j)*x[j];
    }
    for (unsigned j=(i+1);j<n_dof;j++)
    {
     dummy-=(*(matrix_pt))(i,j)*x[j];
    }
    x[i]=dummy/(*(matrix_pt))(i,i);
   } // for(unsigned i=0;i<n_dof;i++)

   // Increment the value of Iterations
   Iterations++;

   // Calculate the residual only if we're not inside the multigrid solver
   if (!Use_as_smoother)
   {
    // Get residual
    matrix_pt->residual(x,rhs,current_residual);

    // Calculate the relative residual norm (i.e. \frac{\|r_{i}\|}{\|r_{0}\|})
    norm_res=current_residual.norm()/norm_f;

    // If required will document convergence history to screen or file (if
    // stream open)
    if (Doc_convergence_history)
    {
     if (!Output_file_stream.is_open())
     {
      oomph_info << Iterations << " " << norm_res << std::endl;
     }
     else
     {
      Output_file_stream << Iterations << " " << norm_res << std::endl;
     }
    } // if (Doc_convergence_history)

    // Check the tolerance only if the residual norm is being computed
    if (norm_res<Tolerance)
    {
     // Break out of the for-loop
     break;
    }
   } // if (!Use_as_smoother)
  } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)

  // Calculate the residual only if we're not inside the multigrid solver
  if (!Use_as_smoother)
  {
   // If time documentation is enabled
   if(Doc_time)
   {
    oomph_info << "\nGS converged. Residual norm: " << norm_res
               << "\nNumber of iterations to convergence: " << Iterations
               << "\n" << std::endl;
   }
  } // if (!Use_as_smoother)

  // Copy result into result
  solution=x;

  // Doc. time for solver
  double t_end = TimingHelpers::timer();
  Solution_time = t_end-t_start;
  if(Doc_time)
  {
   oomph_info << "Time for solve with GS [sec]: "
              << Solution_time << std::endl;
  }

  // If the solver failed to converge and the user asked for an error if
  // this happened
  if((Iterations>Max_iter-1)&&(Throw_error_after_max_iter))
  {
   std::string error_message="Solver failed to converge and you requested ";
   error_message+="an error on convergence failures.";
   throw OomphLibError(error_message,
                       OOMPH_EXCEPTION_LOCATION,
                       OOMPH_CURRENT_FUNCTION);
  }
 } // End of solve_helper


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


//==================================================================
/// \short Explicit template specialisation of the solver for CR
/// matrices: Takes pointer to problem and returns the results
/// vector which contains the solution of the linear system defined
/// by the problem's fully assembled Jacobian and residual vector.
//==================================================================
 void GS<CRDoubleMatrix>::solve(Problem* const &problem_pt,
                                DoubleVector& result)
 {
  // Reset the Use_as_smoother_flag as the solver is not being used
  // as a smoother
  Use_as_smoother=false;
    
  // Find # of degrees of freedom (variables)
  unsigned n_dof=problem_pt->ndof();

  // Initialise timer
  double t_start=TimingHelpers::timer();

  // We're not re-solving
  Resolving=false;

  // Get rid of any previously stored data
  clean_up_memory();

  // Set up the distribution
  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);

  // Build the distribution of the LinearSolver
  this->build_distribution(dist);

  // Allocate space for the Jacobian matrix in CRDoubleMatrix format
  Matrix_pt=new CRDoubleMatrix;

  // Build the matrix
  Matrix_pt->build(this->distribution_pt());

  // Allocate space for the residuals vector
  DoubleVector f;

  // Build it
  f.build(this->distribution_pt(),0.0);

  // Get the global Jacobian and the accompanying residuals vector
  problem_pt->get_jacobian(f,*Matrix_pt);

  // Doc. time for setup
  double t_end=TimingHelpers::timer();
  Jacobian_setup_time=t_end-t_start;

  // Run the generic setup function
  setup_helper(Matrix_pt);

  // We've made the matrix, we can delete it...
  Matrix_can_be_deleted=true;

  // If time documentation is enabled
  if(Doc_time)
  {
   // Output the time for the assembly of the Jacobian to screen
   oomph_info << "Time for setup of Jacobian [sec]: "
              << Jacobian_setup_time << std::endl;
  }

  // Call linear algebra-style solver
  solve_helper(Matrix_pt,f,result);

  // Kill matrix unless it's still required for resolve
  if (!Enable_resolve) clean_up_memory();
 } // End of solve

//==================================================================
/// \short Explicit template specialisation of the smoother_setup
/// function for CR matrices. Set up the smoother for the matrix
/// specified by the pointer. This definition of the smoother_setup
/// has the added functionality that it sorts the entries in the
/// given matrix so that the CRDoubleMatrix implementation of the
/// solver can be used.
//==================================================================
 void GS<CRDoubleMatrix>::setup_helper(DoubleMatrixBase* matrix_pt)
 {
  // Assume the matrix has been passed in from the outside so we must
  // not delete it. This is needed to avoid pre- and post-smoothers
  // deleting the same matrix in the MG solver. If this was originally
  // set to TRUE then this will be sorted out in the other functions
  // from which this was called
  Matrix_can_be_deleted=false;

  // Upcast the input matrix to system matrix to the type MATRIX
  Matrix_pt=dynamic_cast<CRDoubleMatrix*>(matrix_pt);

  // The system matrix here is a CRDoubleMatrix. To make use of the
  // specific implementation of the solver for this type of matrix we
  // need to make sure the entries are arranged correctly
  Matrix_pt->sort_entries();

  // Now get access to the vector Index_of_diagonal_entries
  Index_of_diagonal_entries=Matrix_pt->get_index_of_diagonal_entries();

#ifdef PARANOID
  // Create a boolean variable to store the result of whether or not
  // the entries are in the correct order
  bool is_matrix_sorted=true;

  // Check to make sure the entries have been sorted properly
  is_matrix_sorted=Matrix_pt->entries_are_sorted();

  // If the entries are not sorted properly
  if (!is_matrix_sorted)
  {
   // Throw an error; the solver won't work unless the matrix is sorted
   // properly
   throw OomphLibError("Matrix is not sorted correctly",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif
 } // End of setup_helper

//==================================================================
/// \short Explicit template specialisation of the solve_helper
/// function for CR matrices. Exploiting the sparsity of the given
/// matrix allows for a much faster iterative solver.
//==================================================================
 void GS<CRDoubleMatrix>::solve_helper(DoubleMatrixBase* const &matrix_pt,
                                       const DoubleVector &rhs,
                                       DoubleVector &solution)
 {
  // Get number of dofs
  unsigned n_dof=rhs.nrow();

#ifdef PARANOID
  // Upcast the matrix to the appropriate type
  CRDoubleMatrix* cr_matrix_pt=dynamic_cast<CRDoubleMatrix*>(matrix_pt);
  
  // PARANOID Run the self-tests to check the inputs are correct
  this->check_validity_of_solve_helper_inputs<CRDoubleMatrix>(cr_matrix_pt,rhs,
                                                              solution,n_dof);

  // We don't need the pointer any more but we do still need the matrix
  // so just make tmp_matrix_pt a null pointer
  cr_matrix_pt=0;
#endif

  // Setup the solution if it is not
  if (!solution.distribution_pt()->built())
  {
   solution.build(this->distribution_pt(),0.0);
  }
  // If the solution is already set up
  else
  {
   // If we're inside the multigrid solver then as we traverse up the
   // hierarchy we use the smoother on the updated approximate solution.
   // As such, we should ONLY be resetting all the values to zero if
   // we're NOT inside the multigrid solver
   if (!Use_as_smoother)
   {
    // Initialise the vector with all entries set to zero
    solution.initialise(0.0);
   }
  } // if (!solution.distribution_pt()->built())
  
  // Initialise timer
  double t_start=TimingHelpers::timer();

  // Copy the solution vector into x
  DoubleVector x(solution);
  
  // Create a vector to hold the residual. This will only be built if
  // we're not inside the multigrid solver
  DoubleVector current_residual;

  // Variable to hold the current residual norm. Only used if we're
  // not inside the multigrid solver
  double norm_res=0.0;

  // Variables to hold the initial residual norm. Only used if we're
  // not inside the multigrid solver
  double norm_f=0.0;

  // Initialise the value of Iterations
  Iterations=0;

  // Calculate the residual only if we're not inside the multigrid solver
  if (!Use_as_smoother)
  {
   // Build the residual vector
   current_residual.build(this->distribution_pt(),0.0);
   
   // Calculate the residual r=b-Ax
   matrix_pt->residual(x,rhs,current_residual);

   // Calculate the 2-norm of the residual vector
   norm_res=current_residual.norm();

   // Store the initial norm
   norm_f=norm_res;
   
   // If required will document convergence history to screen
   // or file (if stream is open)
   if (Doc_convergence_history)
   {
    if (!Output_file_stream.is_open())
    {
     oomph_info << Iterations << " " << norm_res << std::endl;
    }
    else
    {
     Output_file_stream << Iterations << " " << norm_res << std::endl;
    }
   } // if (!Use_as_smoother)
  } // if (Doc_convergence_history)

  // In each iteration of Gauss-Seidel we need to calculate
  //              x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
  // where L represents the strict lower triangular portion of A,
  // D represents the diagonal of A and U represents the strict upper
  // triangular portion of A. Since the value of inv(L+D)*rhs remains
  // constant throughout we shall calculate it now and store it to
  // avoid calculating it again with each iteration.

  // Create a temporary matrix pointer to allow for the extraction of the
  // row start, column index and value pointers
  CRDoubleMatrix* tmp_matrix_pt=dynamic_cast<CRDoubleMatrix*>(matrix_pt);
  
  // First acquire access to the value, row_start and column_index arrays
  // from the compressed row matrix
  const double* value_pt=tmp_matrix_pt->value();
  const int* row_start_pt=tmp_matrix_pt->row_start();
  const int* column_index_pt=tmp_matrix_pt->column_index();

  // We've finished using the temporary matrix pointer so make it a null pointer
  tmp_matrix_pt=0;
  
  // We can only calculate this constant term if the diagonal entries
  // of the matrix are nonzero so we check this first
  for (unsigned i=0;i<n_dof;i++)
  {
   // Get the index of the last entry below or on the diagonal in row i
   unsigned diag_index=Index_of_diagonal_entries[i];

   if (unsigned(*(column_index_pt+diag_index))!=i)
   {
    std::string err_strng="Gauss-Seidel cannot operate on a matrix with ";
    err_strng+="zero diagonal entries.";
    throw OomphLibError(err_strng,
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  } // for (unsigned i=0;i<n_dof;i++)

  // If there were no zero diagonal entries we can calculate the vector
  //                 constant_term=inv(L+D)*rhs.

  // Create a vector to hold the constant term in the iteration
  DoubleVector constant_term(this->distribution_pt(),0.0);
  
  // Copy the entries of rhs into the vector, constant_term
  constant_term=rhs;

  // Start by looping over the rows of constant_term
  for (unsigned i=0;i<n_dof;i++)
  {
   // Get the index of the last entry below or on the diagonal in row i
   unsigned diag_index=Index_of_diagonal_entries[i];

   // Get the index of the first entry in row i
   unsigned row_i_start=*(row_start_pt+i);

   // If there are entries strictly below the diagonal then the
   // column index of the first entry in the i-th row cannot be
   // i (since we have ensured that nonzero entries exist on the
   // diagonal so this is the largest the column index of the
   // first entry in this row could be)
   if (unsigned(*(column_index_pt+row_i_start))!=i)
   {
    // Initialise the column index variable
    unsigned column_index=0;

    // Loop over the entries below the diagonal on row i
    for (unsigned j=row_i_start;j<diag_index;j++)
    {
     // Find the column index of this nonzero entry
     column_index=*(column_index_pt+j);

     // Subtract the value of A(i,j)*constant_term(j) (where j<i)
     constant_term[i]-=(*(value_pt+j))*constant_term[column_index];
    }
   } // if (*(column_index_pt+row_i_start)!=i)

   // Finish off by calculating constant_term(i)/A(i,i)
   constant_term[i]/=*(value_pt+diag_index);
  } // for (unsigned i=0;i<n_dof;i++)

  // Build a temporary vector to store the value of A*x with each iteration
  DoubleVector temp_vec(this->distribution_pt(),0.0);

  // Outermost loop: Run max_iter times (the iteration number)
  for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
  {
   // With each iteration we need to calculate
   //            x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
   //                 =-inv(L+D)*U*x_old+constant_term.
   // Since we've already calculated the vector constant_term, it remains
   // for us to calculate inv(L+D)*U*x. To do this we break the calculation
   // into two parts:
   //        (1) The matrix-vector multiplication temp_vec=U*x, and;
   //        (2) The matrix-vector multiplication inv(L+D)*temp_vec.

   //------------------------------
   // (1) Calculating temp_vec=U*x:
   //------------------------------

   // Auxiliary variable to store the index of the first entry on the
   // upper triangular portion of the matrix
   unsigned upper_tri_start=0;

   // Auxiliary variable to store the index of the first entry in the
   // next row
   unsigned next_row_start=0;

   // Set the temporary vector temp_vec to be initially be the zero vector
   temp_vec.initialise(0.0);

   // Loop over the rows of temp_vec (note: we do not loop over the final
   // row since all the entries on the last row of a strict upper triangular
   // matrix are zero)
   for (unsigned i=0;i<n_dof-1;i++)
   {
    // Get the index of the first entry on the upper triangular portion of
    // the matrix noting that the i-th entry in Index_of_diagonal_entries
    // will store the index of the diagonal entry of the i-th row
    upper_tri_start=Index_of_diagonal_entries[i]+1;

    // Get the index of the first entry in the next row
    next_row_start=*(row_start_pt+i+1);

    // Variable to store the column index of each entry
    unsigned column_index=0;

    // Loop over all of the entries above the diagonal
    for (unsigned j=upper_tri_start;j<next_row_start;j++)
    {
     // Get the column index of this entry
     column_index=*(column_index_pt+j);

     // Update the value of temp_vec[i] by adding A(i,j)*x(j) (where j>i)
     temp_vec[i]+=(*(value_pt+j))*x[column_index];
    }
   } // for (unsigned i=0;i<n_dof;i++)

   //-----------------------------------
   // (2) Calculating inv(L+D)*temp_vec:
   //-----------------------------------

   // Now U*x=temp_vec has been calculated we need to calculate the vector
   //                 y=inv(L+D)*U*x=inv(L+D)*temp_vec.
   // Note: the i-th entry in y is given by
   //       y(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*y(j))/A(i,i).
   // We can see from this that to calculate the current entry in y we
   // use its previously calculated entries and the current entry in
   // temp_vec. As a result, we do not need two separate vectors y and
   // temp_vec for this calculation. Instead, we can just update the
   // entries in temp_vec (carefully!) using:
   //   temp_vec(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*temp_vec(j))/A(i,i).

   // Start by looping over the rows of rhs
   for (unsigned i=0;i<n_dof;i++)
   {
    // Get the index of the last entry below or on the diagonal in row i
    unsigned diag_index=Index_of_diagonal_entries[i];

    // Get the index of the first entry in row i
    unsigned row_i_start=*(row_start_pt+i);

    // If there are no entries strictly below the diagonal then the
    // column index of the first entry in the i-th row will be greater
    // than or equal to i
    if (unsigned(*(column_index_pt+row_i_start))!=i)
    {
     // Initialise the column index variable
     unsigned column_index=0;

     // Loop over the entries below the diagonal
     for (unsigned j=row_i_start;j<diag_index;j++)
     {
      // Find the column index of this nonzero entry
      column_index=*(column_index_pt+j);

      // Subtract the value of A(i,j)*y(j) (where j<i)
      temp_vec[i]-=(*(value_pt+j))*temp_vec[column_index];
     }
    } // if (*(column_index_pt+row_i_start)!=i)

    // Finish off by calculating temp_vec(i)/A(i,i)
    temp_vec[i]/=*(value_pt+diag_index);
   } // for (unsigned i=0;i<n_dof;i++)

   //-------------------------
   // Updating the solution x:
   //-------------------------
   
   // The updated solution is given by
   //            x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
   //                 =-temp_vec+constant_term,
   // We have calculated both constant_term and temp_vec so we simply
   // need to calculate x now

   // Set x to be the vector, constant_term
   x=constant_term;

   // Update x to give the next iterate
   x-=temp_vec;

   // Increment the value of Iterations
   Iterations++;
   
   // Calculate the residual only if we're not inside the multigrid solver
   if (!Use_as_smoother)
   {
    // Get residual
    matrix_pt->residual(x,rhs,current_residual);

    // Calculate the relative residual norm (i.e. \frac{\|r_{i}\|}{\|r_{0}\|})
    norm_res=current_residual.norm()/norm_f;

    // If required will document convergence history to screen or file (if
    // stream open)
    if (Doc_convergence_history)
    {
     if (!Output_file_stream.is_open())
     {
      oomph_info << iter_num << " "
                 << norm_res << std::endl;
     }
     else
     {
      Output_file_stream << iter_num << " "
                         << norm_res << std::endl;
     }
    } // if (Doc_convergence_history)
   
    // Check the tolerance only if the residual norm is being computed
    if (norm_res<Tolerance)
    {
     // Break out of the for-loop
     break;
    }
   } // if (!Use_as_smoother)
  } // for (unsigned iter_num=0;iter_num<max_iter;iter_num++)

  // Calculate the residual only if we're not inside the multigrid solver
  if (!Use_as_smoother)
  {
   // If time documentation is enabled
   if(Doc_time)
   {
    oomph_info << "\nGS converged. Residual norm: " << norm_res
               << "\nNumber of iterations to convergence: " << Iterations
               << "\n" << std::endl;
   }
  } // if (!Use_as_smoother)
  
  // Copy the solution into the solution vector
  solution=x;

  // Document the time for the solver
  double t_end=TimingHelpers::timer();
  Solution_time=t_end-t_start;

  // If time documentation is enabled
  if(Doc_time)
  {
   oomph_info << "Time for solve with Gauss-Seidel [sec]: "
              << Solution_time << std::endl;
  }
 
  // If the solver failed to converge and the user asked for an error if
  // this happened
  if((Iterations>Max_iter-1)&&(Throw_error_after_max_iter))
  {
   std::string error_message="Solver failed to converge and you requested ";
   error_message+="an error on convergence failures.";
   throw OomphLibError(error_message,
                       OOMPH_EXCEPTION_LOCATION,
                       OOMPH_CURRENT_FUNCTION);
  }
 } // End of solve_helper


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

 
//==================================================================
/// \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.
//==================================================================
 template<typename MATRIX>
 void DampedJacobi<MATRIX>::solve(Problem* const &problem_pt,
                                  DoubleVector &result)
 {
  // Reset the Use_as_smoother_flag as the solver is not being used
  // as a smoother
  Use_as_smoother=false;
  
  // Find the # of degrees of freedom (variables)
  unsigned n_dof=problem_pt->ndof();

  // Initialise timer
  double t_start=TimingHelpers::timer();

  // We're not re-solving
  Resolving=false;

  // Get rid of any previously stored data
  clean_up_memory();

  // Set up the distribution
  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),
                                 n_dof,false);

  // Assign the distribution to the LinearSolver
  this->build_distribution(dist);

  // Allocate space for the Jacobian matrix in format specified
  // by template parameter
  Matrix_pt=new MATRIX;

  // Get the nonlinear residuals vector
  DoubleVector f;

  // Assign the Jacobian and the residuals vector
  problem_pt->get_jacobian(f,*Matrix_pt);

  // Extract the diagonal entries of the matrix and store them
  extract_diagonal_entries(Matrix_pt);

  // We've made the matrix, we can delete it...
  Matrix_can_be_deleted=true;

  // Doc time for setup
  double t_end=TimingHelpers::timer();
  Jacobian_setup_time=t_end-t_start;

  // If time documentation is enabled
  if(Doc_time)
  {
   oomph_info << "Time for setup of Jacobian [sec]: "
              << Jacobian_setup_time << std::endl;
  }

  // Call linear algebra-style solver
  solve_helper(Matrix_pt,f,result);

  // Kill matrix unless it's still required for resolve
  if (!Enable_resolve) clean_up_memory();
 } // End of solve

//==================================================================
/// \short Linear-algebra-type solver: Takes pointer to a matrix and
/// rhs vector and returns the solution of the linear system.
//==================================================================
 template<typename MATRIX>
 void DampedJacobi<MATRIX>::solve_helper(DoubleMatrixBase* const &matrix_pt,
                                         const DoubleVector& rhs,
                                         DoubleVector& solution)
 {
  // Get number of dofs
  unsigned n_dof=rhs.nrow();

#ifdef PARANOID
  // Upcast the matrix to the appropriate type
  MATRIX* tmp_matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
  
  // PARANOID Run the self-tests to check the inputs are correct
  this->check_validity_of_solve_helper_inputs<MATRIX>(tmp_matrix_pt,rhs,
                                                      solution,n_dof);

  // We don't need the pointer any more but we do still need the matrix
  // so just make tmp_matrix_pt a null pointer
  tmp_matrix_pt=0;
#endif

  // Setup the solution if it is not
  if (!solution.distribution_pt()->built())
  {
   // Build the distribution of the solution vector if it hasn't been done yet
   solution.build(this->distribution_pt(),0.0);
  }
  // If the solution has already been set up
  else
  {   
   // If we're inside the multigrid solver then as we traverse up the
   // hierarchy we use the smoother on the updated approximate solution.
   // As such, we should ONLY be resetting all the values to zero if
   // we're NOT inside the multigrid solver
   if (!Use_as_smoother)
   {
    // Initialise the vector with all entries set to zero
    solution.initialise(0.0);
   }
  } // if (!solution.distribution_pt()->built())

  // Initialise timer
  double t_start = TimingHelpers::timer();

  // Initial guess isn't necessarily zero (restricted solution from finer
  // grids) therefore x needs to be assigned values from the input.
  DoubleVector x(solution);

  // Create a vector to store the value of the constant term, omega*inv(D)*r
  DoubleVector constant_term(this->distribution_pt(),0.0);
 
  // Calculate the constant term vector
  for(unsigned i=0;i<n_dof;i++)
  {
   // Assign the i-th entry of constant_term
   constant_term[i]=Omega*Matrix_diagonal[i]*rhs[i];
  }

  // Create a vector to hold the residual. This will only be built if
  // we're not inside the multigrid solver
  DoubleVector local_residual;

  // Variable to store the 2-norm of the residual vector. Only used
  // if we are not working inside the MG solver
  double norm_res=0.0;

  // Variables to hold the initial residual norm. Only used if we're
  // not inside the multigrid solver
  double norm_f=0.0;
  
  // Initialise the value of Iterations
  Iterations=0;
  
  // Calculate the residual only if we're not inside the multigrid solver
  if (!Use_as_smoother)
  {
   // Build the local residual vector
   local_residual.build(this->distribution_pt(),0.0);

   // Calculate the residual vector
   matrix_pt->residual(x,rhs,local_residual);

   // Calculate the 2-norm
   norm_res=local_residual.norm();

   // Store the initial norm
   norm_f=norm_res;
   
   // If required will document convergence history to screen
   // or file (if stream is open)
   if (Doc_convergence_history)
   {
    if (!Output_file_stream.is_open())
    {
     oomph_info << Iterations << " " << norm_res << std::endl;
    }
    else
    {
     Output_file_stream << Iterations << " " << norm_res <<std::endl;
    }
   } // if (Doc_convergence_history)
  } // if (!Use_as_smoother)

  // Create a temporary vector to store the value of A*e
  // on each iteration and another to store the residual
  // at the end of each iteration
  DoubleVector temp_vec(this->distribution_pt(),0.0);

  // Outermost loop: Run up to Max_iter times (the iteration number)
  for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
  {
   // Calculate A*e
   matrix_pt->multiply(x,temp_vec);

   // Loop over each degree of freedom and update
   // the current approximation
   for (unsigned idof=0;idof<n_dof;idof++)
   {
    // Scale the idof'th entry of temp_vec
    // by omega/A(i,i)
    temp_vec[idof]*=Omega*Matrix_diagonal[idof];
   }

   // Update the value of e (in the system Ae=r)
   x+=constant_term;
   x-=temp_vec;

   // Increment the value of Iterations
   Iterations++;
   
   // Calculate the residual only if we're not inside the multigrid solver
   if (!Use_as_smoother)
   {
    // Get residual
    matrix_pt->residual(x,rhs,local_residual);

    // Calculate the 2-norm of the residual r=b-Ax
    norm_res=local_residual.norm()/norm_f;

    // If required, this will document convergence history to
    // screen or file (if the stream is open)
    if (Doc_convergence_history)
    {
     if (!Output_file_stream.is_open())
     {
      oomph_info << Iterations << " "
                 << norm_res << std::endl;
     }
     else
     {
      Output_file_stream << Iterations << " "
                         << norm_res << std::endl;
     }
    } // if (Doc_convergence_history)

    // Check the tolerance only if the residual norm is being computed
    if (norm_res<Tolerance)
    {
     // Break out of the for-loop
     break;
    }
   } // if (!Use_as_smoother)
  } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)

  // Calculate the residual only if we're not inside the multigrid solver
  if (!Use_as_smoother)
  {
   // If time documentation is enabled
   if(Doc_time)
   {
    oomph_info << "\nDamped Jacobi converged. Residual norm: " << norm_res
               << "\nNumber of iterations to convergence: " << Iterations
               << "\n" << std::endl;
   }
  } // if (!Use_as_smoother)
  
  // Copy the solution into the solution vector
  solution=x;

  // Doc. time for solver
  double t_end=TimingHelpers::timer();
  Solution_time=t_end-t_start;
  if(Doc_time)
  {
   oomph_info << "Time for solve with damped Jacobi [sec]: "
              << Solution_time << std::endl;
  }
 
  // If the solver failed to converge and the user asked for an error if
  // this happened
  if((Iterations>Max_iter-1)&&(Throw_error_after_max_iter))
  {
   std::string error_message="Solver failed to converge and you requested ";
   error_message+="an error on convergence failures.";
   throw OomphLibError(error_message,
                       OOMPH_EXCEPTION_LOCATION,
                       OOMPH_CURRENT_FUNCTION);
  }
 } // End of solve_helper


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


//==================================================================
/// \Short Re-solve the system defined by the last assembled Jacobian
/// and the rhs vector specified here. Solution is returned in
/// the vector result.
//==================================================================
 template<typename MATRIX>
 void GMRES<MATRIX>::resolve(const DoubleVector &rhs,
                             DoubleVector &result)
 {
  // We are re-solving
  Resolving=true;

#ifdef PARANOID
  if (Matrix_pt==0)
  {
   throw OomphLibError("No matrix was stored -- cannot re-solve",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

  // Call linear algebra-style solver
  this->solve(Matrix_pt,rhs,result);

  // Reset re-solving flag
  Resolving=false;
 }


//==================================================================
/// Solver: Takes pointer to problem and returns the results vector
/// which contains the solution of the linear system defined by
/// the problem's fully assembled Jacobian and residual vector.
//==================================================================
 template<typename MATRIX>
 void GMRES<MATRIX>::solve(Problem* const &problem_pt, DoubleVector &result)
 {

  //Find # of degrees of freedom (variables)
  unsigned n_dof = problem_pt->ndof();

  // Initialise timer
  double t_start = TimingHelpers::timer();

  // We're not re-solving
  Resolving=false;

  // Get rid of any previously stored data
  clean_up_memory();

  // setup the distribution
  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),
                                 n_dof,false);
  this->build_distribution(dist);

  // Get Jacobian matrix in format specified by template parameter
  // and nonlinear residual vector
  Matrix_pt=new MATRIX;
  DoubleVector f;
  if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt) != 0)
  {
   if (dynamic_cast<CRDoubleMatrix*>(Matrix_pt) != 0)
   {
    dynamic_cast<CRDoubleMatrix* >(Matrix_pt)->build(this->distribution_pt());
    f.build(this->distribution_pt(),0.0);
   }
  }
  problem_pt->get_jacobian(f,*Matrix_pt);

  // We've made the matrix, we can delete it...
  Matrix_can_be_deleted=true;

  // Doc time for setup
  double t_end = TimingHelpers::timer();
  Jacobian_setup_time= t_end-t_start;

  if(Doc_time)
  {
   oomph_info << "Time for setup of Jacobian [sec]: "
              << Jacobian_setup_time << std::endl;
  }

  // Call linear algebra-style solver
  //If the result distribution is wrong, then redistribute
  //before the solve and return to original distribution
  //afterwards
  if((result.built()) &&
     (!(*result.distribution_pt() == *this->distribution_pt())))
  {
   LinearAlgebraDistribution
    temp_global_dist(result.distribution_pt());
   result.build(this->distribution_pt(),0.0);
   this->solve_helper(Matrix_pt,f,result);
   result.redistribute(&temp_global_dist);
  }
  //Otherwise just solve
  else
  {
   this->solve_helper(Matrix_pt,f,result);
  }

  // Kill matrix unless it's still required for resolve
  if (!Enable_resolve) clean_up_memory();

 };


//=============================================================================
/// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
/// and returns the solution of the linear system.
/// based on the algorithm presented in Templates for the
/// Solution of Linear Systems: Building Blocks for Iterative Methods, Barrett,
/// Berry et al, SIAM, 2006 and the implementation in the IML++ library :
/// http://math.nist.gov/iml++/
//=============================================================================
 template <typename MATRIX>
 void GMRES<MATRIX>::solve_helper(DoubleMatrixBase* const &matrix_pt,
                                  const DoubleVector &rhs,
                                  DoubleVector &solution)
 {

  // Get number of dofs
  unsigned n_dof=rhs.nrow();

#ifdef PARANOID
  // PARANOID check that if the matrix is distributable then it should not be
  // then it should not be distributed
  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
  {
   if (dynamic_cast<DistributableLinearAlgebraObject*>
       (matrix_pt)->distributed())
   {
    std::ostringstream error_message_stream;
    error_message_stream
     << "The matrix must not be distributed.";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  }
  // PARANOID check that this rhs distribution is setup
  if (!rhs.built())
  {
   std::ostringstream error_message_stream;
   error_message_stream
    << "The rhs vector distribution must be setup.";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
  // PARANOID check that the rhs has the right number of global rows
  if(rhs.nrow() != n_dof)
  {
   throw OomphLibError(
    "RHS does not have the same dimension as the linear system",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
  // PARANOID check that the rhs is not distributed
  if (rhs.distribution_pt()->distributed())
  {
   std::ostringstream error_message_stream;
   error_message_stream
    << "The rhs vector must not be distributed.";
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
  // PARANOID check that if the result is setup it matches the distribution
  // of the rhs
  if (solution.built())
  {
   if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
   {
    std::ostringstream error_message_stream;
    error_message_stream
     << "If the result distribution is setup then it must be the same as the "
     << "rhs distribution";
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  }
#endif

  // Set up the solution if it is not
  if (!solution.built())
  {
   solution.build(this->distribution_pt(),0.0);
  }
  // Otherwise initialise to zero
  else
  {
   solution.initialise(0.0);
  }

  // Time solver
  double t_start=TimingHelpers::timer();

  // Relative residual
  double resid;

  // iteration counter
  unsigned iter = 1;

  // if not using iteration restart set Restart to n_dof
  if (!Iteration_restart)
  {
   Restart = n_dof;
  }

  // initialise vectors
  Vector<double> s(Restart + 1,0);
  Vector<double> cs(Restart + 1);
  Vector<double> sn(Restart + 1);
  DoubleVector w(this->distribution_pt(),0.0);

  // Setup preconditioner only if we're not re-solving
  if (!Resolving)
  {
   // only setup the preconditioner before solve if require
   if (Setup_preconditioner_before_solve)
   {

    //Setup preconditioner from the Jacobian matrix
    double t_start_prec = TimingHelpers::timer();

    // do not setup
    preconditioner_pt()->setup(matrix_pt);

    // Doc time for setup of preconditioner
    double t_end_prec = TimingHelpers::timer();
    Preconditioner_setup_time = t_end_prec-t_start_prec;

    if(Doc_time)
    {
     oomph_info << "Time for setup of preconditioner  [sec]: "
                << Preconditioner_setup_time << std::endl;
    }
   }
  }
  else
  {
   if(Doc_time)
   {
    oomph_info << "Setup of preconditioner is bypassed in resolve mode"
               << std::endl;
   }
  }

  // solve b-Jx = Mr for r (assumes x = 0);
  DoubleVector r(this->distribution_pt(),0.0);
  if(Preconditioner_LHS)
  {
   preconditioner_pt()->preconditioner_solve(rhs,r);
  }
  else
  {
   r=rhs;
  }
  double normb = 0;

  // compute norm(r)
  double* r_pt = r.values_pt();
  for (unsigned i = 0; i < n_dof; i++)
  {
   normb += r_pt[i] * r_pt[i];
  }
  normb = sqrt(normb);

  // set beta (the initial residual)
  double beta = normb;

  // compute initial relative residual
  if (normb == 0.0) normb = 1;
  resid = beta / normb;

  // if required will document convergence history to screen or file (if
  // stream open)
  if (Doc_convergence_history)
  {
   if (!Output_file_stream.is_open())
   {
    oomph_info <<  0 << " " << resid << std::endl;
   }
   else
   {
    Output_file_stream << 0 << " " << resid <<std::endl;
   }
  }

  // if GMRES converges immediately
  if (resid <= Tolerance)
  {
   if(Doc_time)
   {
    oomph_info << "GMRES converged immediately. Normalised residual norm: "
               << resid << std::endl;
   }
   // Doc time for solver
   double t_end = TimingHelpers::timer();
   Solution_time = t_end-t_start;

   if(Doc_time)
   {
    oomph_info << "Time for solve with GMRES  [sec]: "
               << Solution_time << std::endl;
   }
   return;
  }


  // initialise vector of orthogonal basis vectors (v) and upper hessenberg
  // matrix H
  // NOTE: for implementation purpose the upper hessenberg matrix indexes are
  // are swapped so the matrix is effectively transposed
  Vector<DoubleVector> v;
  v.resize(Restart + 1);
  Vector<Vector<double> > H(Restart + 1);

  // while...
  while (iter <= Max_iter)
  {

   // set zeroth basis vector v[0] to r/beta
   v[0].build(this->distribution_pt(),0.0);
   double* v0_pt = v[0].values_pt();
   for (unsigned i = 0; i < n_dof; i++)
   {
    v0_pt[i] = r_pt[i] / beta;
   }

   //
   s[0] = beta;

   // inner iteration counter for restarted version
   unsigned iter_restart;

   // perform iterations
   for (iter_restart = 0; iter_restart < Restart && iter <= Max_iter;
        iter_restart++, iter++)
   {

    // resize next column of upper hessenberg matrix
    H[iter_restart].resize(iter_restart+2);

    // solve Jv[i] = Mw for w
    {
     DoubleVector temp(this->distribution_pt(),0.0);
     if(Preconditioner_LHS)
     {
      // solve Jv[i] = Mw for w
      matrix_pt->multiply(v[iter_restart],temp);
      preconditioner_pt()->preconditioner_solve(temp,w);
     }
     else
     {
      //w=JM^{-1}v by saad p270
      preconditioner_pt()->preconditioner_solve(v[iter_restart],temp);
      matrix_pt->multiply(temp,w);
     }
    }

    //
    double* w_pt = w.values_pt();
    for (unsigned k = 0; k <= iter_restart; k++)
    {

     //
     H[iter_restart][k] = 0.0;
     double* vk_pt = v[k].values_pt();
     for (unsigned i = 0; i < n_dof; i++)
     {
      H[iter_restart][k] += w_pt[i]*vk_pt[i];
     }

     //
     for (unsigned i = 0; i < n_dof; i++)
     {
      w_pt[i] -= H[iter_restart][k] * vk_pt[i];
     }
    }

    //
    {
     double temp_norm_w = 0.0;
     for (unsigned i = 0; i < n_dof; i++)
     {
      temp_norm_w += w_pt[i]*w_pt[i];
     }
     temp_norm_w = sqrt(temp_norm_w);
     H[iter_restart][iter_restart+1] = temp_norm_w;
    }

    //
    v[iter_restart + 1].build(this->distribution_pt(),0.0);
    double* v_pt = v[iter_restart + 1].values_pt();
    for (unsigned i = 0; i < n_dof; i++)
    {
     v_pt[i] = w_pt[i] / H[iter_restart][iter_restart+1];
    }

    //
    for (unsigned k = 0; k < iter_restart; k++)
    {
     apply_plane_rotation(H[iter_restart][k], H[iter_restart][k+1], cs[k],
                          sn[k]);
    }
    generate_plane_rotation(H[iter_restart][iter_restart],
                            H[iter_restart][iter_restart+1],
                            cs[iter_restart],
                            sn[iter_restart]);
    apply_plane_rotation(H[iter_restart][iter_restart],
                         H[iter_restart][iter_restart+1], cs[iter_restart],
                         sn[iter_restart]);
    apply_plane_rotation(s[iter_restart],s[iter_restart+1],cs[iter_restart],
                         sn[iter_restart]);

    // compute current residual
    beta = std::fabs(s[iter_restart+1]);

    // compute relative residual
    resid = beta/normb;

    // if required will document convergence history to screen or file (if
    // stream open)
    if (Doc_convergence_history)
    {
     if (!Output_file_stream.is_open())
     {
      oomph_info << iter << " "<< resid << std::endl;
     }
     else
     {
      Output_file_stream << iter << " " << resid <<std::endl;
     }
    }

    // if required tolerance found
    if (resid < Tolerance)
    {
     // update result vector
     update(iter_restart, H, s, v, solution);

     // document convergence
     if(Doc_time)
     {
      oomph_info << std::endl;
      oomph_info << "GMRES converged (1). Normalised residual norm: "
                 << resid << std::endl;
      oomph_info << "Number of iterations to convergence: "
                 << iter << std::endl;
      oomph_info << std::endl;
     }

     // Doc time for solver
     double t_end = TimingHelpers::timer();
     Solution_time = t_end-t_start;

     Iterations = iter;

     if(Doc_time)
     {
      oomph_info << "Time for solve with GMRES  [sec]: "
                 << Solution_time << std::endl;
     }
     return;
    }
   }

   // update
   if (iter_restart>0) update((iter_restart - 1), H, s, v, solution);

   // solve Mr = (b-Jx) for r
   {
    DoubleVector temp(this->distribution_pt(),0.0);
    matrix_pt->multiply(solution,temp);
    double* temp_pt = temp.values_pt();
    const double* rhs_pt = rhs.values_pt();
    for (unsigned i = 0; i < n_dof; i++)
    {
     temp_pt[i] = rhs_pt[i] - temp_pt[i];
    }

    if(Preconditioner_LHS)
    {
     preconditioner_pt()->preconditioner_solve(temp,r);
    }
   }

   // compute current residual
   beta = 0.0;
   r_pt = r.values_pt();
   for (unsigned i = 0; i < n_dof; i++)
   {
    beta += r_pt[i] * r_pt[i];
   }
   beta = sqrt(beta);

   // if relative residual within tolerance
   resid = beta / normb;
   if (resid < Tolerance)
   {
    if(Doc_time)
    {
     oomph_info << std::endl;
     oomph_info << "GMRES converged (2). Normalised residual norm: "
                << resid << std::endl;
     oomph_info << "Number of iterations to convergence: "
                << iter << std::endl;
     oomph_info << std::endl;
    }

    // Doc time for solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end-t_start;

    Iterations = iter;

    if(Doc_time)
    {
     oomph_info << "Time for solve with GMRES  [sec]: "
                << Solution_time << std::endl;
    }
    return;
   }
  }


  // otherwise GMRES failed convergence
  oomph_info << std::endl;
  oomph_info << "GMRES did not converge to required tolerance! "
             << std::endl;
  oomph_info << "Returning with normalised residual norm: "
             << resid << std::endl;
  oomph_info << "after " << Max_iter<< " iterations." << std::endl;
  oomph_info << std::endl;


  if(Throw_error_after_max_iter)
  {
   std::string err = "Solver failed to converge and you requested an error";
   err += " on convergence failures.";
   throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                       OOMPH_CURRENT_FUNCTION);
  }

  return;

 } // End GMRES




//Ensure build of required objects
 template class BiCGStab<CCDoubleMatrix>;
 template class BiCGStab<CRDoubleMatrix>;
 template class BiCGStab<DenseDoubleMatrix>;

 template class CG<CCDoubleMatrix>;
 template class CG<CRDoubleMatrix>;
 template class CG<DenseDoubleMatrix>;

 template class GS<CCDoubleMatrix>;
 template class GS<CRDoubleMatrix>;
 template class GS<DenseDoubleMatrix>;
 
 template class DampedJacobi<CCDoubleMatrix>;
 template class DampedJacobi<CRDoubleMatrix>;
 template class DampedJacobi<DenseDoubleMatrix>;

 template class GMRES<CCDoubleMatrix>;
 template class GMRES<CRDoubleMatrix>;
 template class GMRES<DenseDoubleMatrix>;

// Solvers for SumOfMatrices class
 template class BiCGStab<SumOfMatrices>;
 template class CG<SumOfMatrices>;
 template class GS<SumOfMatrices>;
 template class GMRES<SumOfMatrices>;


}