The GMRES method. More...

#include <iterative_linear_solver.h>

Inheritance diagram for oomph::GMRES< MATRIX >:

Public Member Functions
	GMRES ()
	Constructor. More...

virtual	~GMRES ()
	Destructor (cleanup storage) More...

	GMRES (const GMRES &)=delete
	Broken copy constructor. More...

void	operator= (const GMRES &)=delete
	Broken assignment operator. More...

void	disable_resolve ()
	Overload disable resolve so that it cleans up memory too. More...

void	enable_computation_of_gradient ()
	function to enable the computation of the gradient More...

void	solve (Problem *const &problem_pt, DoubleVector &result)

void	solve (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)

void	solve (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)

void	resolve (const DoubleVector &rhs, DoubleVector &result)

unsigned	iterations () const
	Number of iterations taken. More...

bool	iteration_restart () const
	access function indicating whether restarted GMRES is used More...

void	enable_iteration_restart (const unsigned &restart)

void	disable_iteration_restart ()
	switches off iteration restart More...

void	set_preconditioner_LHS ()
	Set left preconditioning (the default) More...

void	set_preconditioner_RHS ()
	Enable right preconditioning. More...

Public Member Functions inherited from oomph::IterativeLinearSolver
	IterativeLinearSolver ()

	IterativeLinearSolver (const IterativeLinearSolver &)=delete
	Broken copy constructor. More...

void	operator= (const IterativeLinearSolver &)=delete
	Broken assignment operator. More...

virtual	~IterativeLinearSolver ()
	Destructor (empty) More...

Preconditioner *&	preconditioner_pt ()
	Access function to preconditioner. More...

Preconditioner *const &	preconditioner_pt () const
	Access function to preconditioner (const version) More...

double &	tolerance ()
	Access to convergence tolerance. More...

unsigned &	max_iter ()
	Access to max. number of iterations. More...

void	enable_doc_convergence_history ()
	Enable documentation of the convergence history. More...

void	disable_doc_convergence_history ()
	Disable documentation of the convergence history. More...

void	open_convergence_history_file_stream (const std::string &file_name, const std::string &zone_title="")

void	close_convergence_history_file_stream ()
	Close convergence history output stream. More...

double	jacobian_setup_time () const

double	linear_solver_solution_time () const
	return the time taken to solve the linear system More...

virtual double	preconditioner_setup_time () const
	returns the the time taken to setup the preconditioner More...

void	enable_setup_preconditioner_before_solve ()
	Setup the preconditioner before the solve. More...

void	disable_setup_preconditioner_before_solve ()
	Don't set up the preconditioner before the solve. More...

void	enable_error_after_max_iter ()
	Throw an error if we don't converge within max_iter. More...

void	disable_error_after_max_iter ()
	Don't throw an error if we don't converge within max_iter (default). More...

void	enable_iterative_solver_as_preconditioner ()

void	disable_iterative_solver_as_preconditioner ()

Public Member Functions inherited from oomph::LinearSolver
	LinearSolver ()
	Empty constructor, initialise the member data. More...

	LinearSolver (const LinearSolver &dummy)=delete
	Broken copy constructor. More...

void	operator= (const LinearSolver &)=delete
	Broken assignment operator. More...

virtual	~LinearSolver ()
	Empty virtual destructor. More...

void	enable_doc_time ()
	Enable documentation of solve times. More...

void	disable_doc_time ()
	Disable documentation of solve times. More...

bool	is_doc_time_enabled () const
	Is documentation of solve times enabled? More...

bool	is_resolve_enabled () const
	Boolean flag indicating if resolves are enabled. More...

virtual void	enable_resolve ()

virtual void	solve_transpose (Problem *const &problem_pt, DoubleVector &result)

virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)

virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)

virtual void	resolve_transpose (const DoubleVector &rhs, DoubleVector &result)

void	disable_computation_of_gradient ()

void	reset_gradient ()

void	get_gradient (DoubleVector &gradient)
	function to access the gradient, provided it has been computed More...

Public Member Functions inherited from oomph::DistributableLinearAlgebraObject
	DistributableLinearAlgebraObject ()
	Default constructor - create a distribution. More...

	DistributableLinearAlgebraObject (const DistributableLinearAlgebraObject &matrix)=delete
	Broken copy constructor. More...

void	operator= (const DistributableLinearAlgebraObject &)=delete
	Broken assignment operator. More...

virtual	~DistributableLinearAlgebraObject ()
	Destructor. More...

LinearAlgebraDistribution *	distribution_pt () const
	access to the LinearAlgebraDistribution More...

unsigned	nrow () const
	access function to the number of global rows. More...

unsigned	nrow_local () const
	access function for the num of local rows on this processor. More...

unsigned	nrow_local (const unsigned &p) const
	access function for the num of local rows on this processor. More...

unsigned	first_row () const
	access function for the first row on this processor More...

unsigned	first_row (const unsigned &p) const
	access function for the first row on this processor More...

bool	distributed () const
	distribution is serial or distributed More...

bool	distribution_built () const

void	build_distribution (const LinearAlgebraDistribution *const dist_pt)

void	build_distribution (const LinearAlgebraDistribution &dist)

Private Member Functions
void	solve_helper (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
	General interface to solve function. More...

void	clean_up_memory ()
	Cleanup data that's stored for resolve (if any has been stored) More...

void	update (const unsigned &k, const Vector< Vector< double >> &H, const Vector< double > &s, const Vector< DoubleVector > &v, DoubleVector &x)

void	generate_plane_rotation (double &dx, double &dy, double &cs, double &sn)

void	apply_plane_rotation (double &dx, double &dy, double &cs, double &sn)

Private Attributes
unsigned	Iterations
	Number of iterations taken. More...

unsigned	Restart

bool	Iteration_restart
	boolean indicating if iteration restarting is used More...

MATRIX *	Matrix_pt
	Pointer to matrix. More...

bool	Resolving

bool	Matrix_can_be_deleted

bool	Preconditioner_LHS

double	Preconditioner_application_time
	Storage for the time spent applying the preconditioner. More...

Additional Inherited Members
Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void	clear_distribution ()

Protected Attributes inherited from oomph::IterativeLinearSolver
bool	Doc_convergence_history

std::ofstream	Output_file_stream
	Output file stream for convergence history. More...

double	Tolerance
	Convergence tolerance. More...

unsigned	Max_iter
	Maximum number of iterations. More...

Preconditioner *	Preconditioner_pt
	Pointer to the preconditioner. More...

double	Jacobian_setup_time
	Jacobian setup time. More...

double	Solution_time
	linear solver solution time More...

double	Preconditioner_setup_time
	Preconditioner setup time. More...

bool	Setup_preconditioner_before_solve

bool	Throw_error_after_max_iter

bool	Use_iterative_solver_as_preconditioner
	Use the iterative solver as preconditioner. More...

bool	First_time_solve_when_used_as_preconditioner

Protected Attributes inherited from oomph::LinearSolver
bool	Enable_resolve

bool	Doc_time
	Boolean flag that indicates whether the time taken. More...

bool	Compute_gradient

bool	Gradient_has_been_computed
	flag that indicates whether the gradient was computed or not More...

DoubleVector	Gradient_for_glob_conv_newton_solve

Static Protected Attributes inherited from oomph::IterativeLinearSolver
static IdentityPreconditioner	Default_preconditioner

Detailed Description

template<typename MATRIX>
class oomph::GMRES< MATRIX >

The GMRES method.

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

Constructor & Destructor Documentation

◆ GMRES() [1/2]

template<typename MATRIX >

oomph::GMRES< MATRIX >::GMRES ( )

inline

Constructor.

       : Iterations(0),
         Matrix_pt(0),
         Resolving(false),
         Matrix_can_be_deleted(true),
         Preconditioner_application_time(0.0)
     {
       Preconditioner_LHS = true;
       Iteration_restart = false;
     }

References oomph::GMRES< MATRIX >::Iteration_restart, and oomph::GMRES< MATRIX >::Preconditioner_LHS.

◆ ~GMRES()

template<typename MATRIX >

virtual oomph::GMRES< MATRIX >::~GMRES ( )

inlinevirtual

Destructor (cleanup storage)

     {
       clean_up_memory();
     }

References oomph::GMRES< MATRIX >::clean_up_memory().

◆ GMRES() [2/2]

template<typename MATRIX >

oomph::GMRES< MATRIX >::GMRES ( const GMRES< MATRIX > & )

delete

Broken copy constructor.

Member Function Documentation

◆ apply_plane_rotation()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::apply_plane_rotation	(	double &	dx,
		double &	dy,
		double &	cs,
		double &	sn
	)

inlineprivate

Helper function: Apply plane rotation. This is done using the update:

\[ \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}. \]

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

◆ clean_up_memory()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::clean_up_memory ( )

inlineprivatevirtual

Cleanup data that's stored for resolve (if any has been stored)

Reimplemented from oomph::LinearSolver.

     {
       if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
       {
         delete Matrix_pt;
         Matrix_pt = 0;
       }
     }

References oomph::GMRES< MATRIX >::Matrix_can_be_deleted, and oomph::GMRES< MATRIX >::Matrix_pt.

Referenced by oomph::GMRES< MATRIX >::disable_resolve(), and oomph::GMRES< MATRIX >::~GMRES().

◆ disable_iteration_restart()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::disable_iteration_restart ( )

inline

switches off iteration restart

     {
       Iteration_restart = false;
     }

References oomph::GMRES< MATRIX >::Iteration_restart.

◆ disable_resolve()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::disable_resolve ( )

inlinevirtual

Overload disable resolve so that it cleans up memory too.

Reimplemented from oomph::LinearSolver.

     {
       LinearSolver::disable_resolve();
       clean_up_memory();
     }

References oomph::GMRES< MATRIX >::clean_up_memory(), and oomph::LinearSolver::disable_resolve().

◆ enable_computation_of_gradient()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::enable_computation_of_gradient ( )

inlinevirtual

function to enable the computation of the gradient

Reimplemented from oomph::LinearSolver.

     {
       Compute_gradient = true;
     }

References oomph::LinearSolver::Compute_gradient.

◆ enable_iteration_restart()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::enable_iteration_restart ( const unsigned & restart )

inline

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

     {
       Restart = restart;
       Iteration_restart = true;
     }

References oomph::GMRES< MATRIX >::Iteration_restart.

◆ generate_plane_rotation()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::generate_plane_rotation	(	double &	dx,
		double &	dy,
		double &	cs,
		double &	sn
	)

inlineprivate

Helper function: Generate a plane rotation. This is done by finding the values of \( \cos(\theta) \) (i.e. cs) and \sin(\theta) (i.e. sn) such that:

\[ \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}, \]

where \( r=\sqrt{pow(dx,2)+pow(dy,2)} \). The values of a and b are given by:

\[ \cos\theta&=\dfrac{dx}{\sqrt{pow(dx,2)+pow(dy,2)}}, \]

and

\[ \sin\theta&=\dfrac{dy}{\sqrt{pow(dx,2)+pow(dy,2)}}. \]

Taken from: Saad Y."Iterative methods for sparse linear systems", p.192

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

References boost::multiprecision::fabs(), and sqrt().

◆ iteration_restart()

template<typename MATRIX >

bool oomph::GMRES< MATRIX >::iteration_restart ( ) const

inline

access function indicating whether restarted GMRES is used

     {
       return Iteration_restart;
     }

References oomph::GMRES< MATRIX >::Iteration_restart.

◆ iterations()

template<typename MATRIX >

unsigned oomph::GMRES< MATRIX >::iterations ( ) const

inlinevirtual

Number of iterations taken.

Implements oomph::IterativeLinearSolver.

     {
       return Iterations;
     }

References oomph::GMRES< MATRIX >::Iterations.

◆ operator=()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::operator= ( const GMRES< MATRIX > & )

delete

Broken assignment operator.

◆ resolve()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::resolve	(	const DoubleVector &	rhs,
		DoubleVector &	result
	)

virtual

Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here. Solution is returned in the vector 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.

Reimplemented from oomph::LinearSolver.

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

References OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and solve.

◆ set_preconditioner_LHS()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::set_preconditioner_LHS ( )

inline

Set left preconditioning (the default)

     {
       Preconditioner_LHS = true;
     }

References oomph::GMRES< MATRIX >::Preconditioner_LHS.

◆ set_preconditioner_RHS()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::set_preconditioner_RHS ( )

inline

Enable right preconditioning.

     {
       Preconditioner_LHS = false;
     }

References oomph::GMRES< MATRIX >::Preconditioner_LHS.

◆ solve() [1/3]

template<typename MATRIX >

void oomph::GMRES< MATRIX >::solve	(	DoubleMatrixBase *const &	matrix_pt,
		const DoubleVector &	rhs,
		DoubleVector &	solution
	)

inlinevirtual

Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the linear system.

Reimplemented from oomph::LinearSolver.

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

References oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Enable_resolve, oomph::GMRES< MATRIX >::Matrix_can_be_deleted, oomph::GMRES< MATRIX >::Matrix_pt, oomph::GMRES< MATRIX >::Resolving, BiharmonicTestFunctions1::solution(), and oomph::GMRES< MATRIX >::solve_helper().

◆ solve() [2/3]

template<typename MATRIX >

void oomph::GMRES< MATRIX >::solve	(	DoubleMatrixBase *const &	matrix_pt,
		const Vector< double > &	rhs,
		Vector< double > &	result
	)

inlinevirtual

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

Reimplemented from oomph::LinearSolver.

     {
       LinearSolver::solve(matrix_pt, rhs, result);
     }

References oomph::LinearSolver::solve().

◆ solve() [3/3]

template<typename MATRIX >

void oomph::GMRES< MATRIX >::solve	(	Problem *const &	problem_pt,
		DoubleVector &	result
	)

virtual

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.

Implements oomph::LinearSolver.

   {
     // Find # of degrees of freedom (variables)
     unsigned n_dof = problem_pt->ndof();
  
     // Initialise timer
     double t_start = TimingHelpers::timer();
  
     // We're not re-solving
     Resolving = false;
  
     // Get rid of any previously stored data
     clean_up_memory();
  
     // setup the distribution
     LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
     this->build_distribution(dist);
  
     // Get Jacobian matrix in format specified by template parameter
     // and nonlinear residual vector
     Matrix_pt = new MATRIX;
     DoubleVector f;
     if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt) != 0)
     {
       if (dynamic_cast<CRDoubleMatrix*>(Matrix_pt) != 0)
       {
         dynamic_cast<CRDoubleMatrix*>(Matrix_pt)->build(
           this->distribution_pt());
         f.build(this->distribution_pt(), 0.0);
       }
     }
     problem_pt->get_jacobian(f, *Matrix_pt);
  
     // We've made the matrix, we can delete it...
     Matrix_can_be_deleted = true;
  
     // Doc time for setup
     double t_end = TimingHelpers::timer();
     Jacobian_setup_time = t_end - t_start;
  
     if (Doc_time)
     {
       oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                  << std::endl;
     }
  
     // If we want to compute the gradient for the globally convergent
     // Newton method, then do it here
     if (Compute_gradient)
     {
       // Compute it
       Matrix_pt->multiply_transpose(f, Gradient_for_glob_conv_newton_solve);
       // Set the flag
       Gradient_has_been_computed = true;
     }
  
     // Call linear algebra-style solver
     // If the result distribution is wrong, then redistribute
     // before the solve and return to original distribution
     // afterwards
     if ((result.built()) &&
         (!(*result.distribution_pt() == *this->distribution_pt())))
     {
       LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
       result.build(this->distribution_pt(), 0.0);
       this->solve_helper(Matrix_pt, f, result);
       result.redistribute(&temp_global_dist);
     }
     // Otherwise just solve
     else
     {
       this->solve_helper(Matrix_pt, f, result);
     }
  
     // Kill matrix unless it's still required for resolve
     if (!Enable_resolve) clean_up_memory();
   };

References oomph::DoubleVector::build(), oomph::DoubleVector::built(), oomph::TerminateHelper::clean_up_memory(), oomph::Problem::communicator_pt(), oomph::DistributableLinearAlgebraObject::distribution_pt(), f(), oomph::Problem::get_jacobian(), oomph::Problem::ndof(), oomph::oomph_info, oomph::DoubleVector::redistribute(), and oomph::TimingHelpers::timer().

◆ solve_helper()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::solve_helper	(	DoubleMatrixBase *const &	matrix_pt,
		const DoubleVector &	rhs,
		DoubleVector &	solution
	)

private

General interface to solve function.

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

   {
     // Get number of dofs
     unsigned n_dof = rhs.nrow();
  
 #ifdef PARANOID
     // PARANOID check that if the matrix is distributable then it should not be
     // then it should not be distributed
     if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
     {
       if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
             ->distributed())
       {
         std::ostringstream error_message_stream;
         error_message_stream << "The matrix must not be distributed.";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
     // PARANOID check that this rhs distribution is setup
     if (!rhs.built())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The rhs vector distribution must be setup.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // PARANOID check that the rhs has the right number of global rows
     if (rhs.nrow() != n_dof)
     {
       throw OomphLibError(
         "RHS does not have the same dimension as the linear system",
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
     }
     // PARANOID check that the rhs is not distributed
     if (rhs.distribution_pt()->distributed())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The rhs vector must not be distributed.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // PARANOID check that if the result is setup it matches the distribution
     // of the rhs
     if (solution.built())
     {
       if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
       {
         std::ostringstream error_message_stream;
         error_message_stream << "If the result distribution is setup then it "
                                 "must be the same as the "
                              << "rhs distribution";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // Reset the time spent applying the preconditioner
     Preconditioner_application_time = 0.0;
  
     // Set up the solution if it is not
     if (!solution.built())
     {
       solution.build(this->distribution_pt(), 0.0);
     }
     // Otherwise initialise to zero
     else
     {
       solution.initialise(0.0);
     }
  
     // Time solver
     double t_start = TimingHelpers::timer();
  
     // Relative residual
     double resid;
  
     // iteration counter
     unsigned iter = 1;
  
     // if not using iteration restart set Restart to n_dof
     if (!Iteration_restart)
     {
       Restart = n_dof;
     }
  
     // initialise vectors
     Vector<double> s(Restart + 1, 0);
     Vector<double> cs(Restart + 1);
     Vector<double> sn(Restart + 1);
     DoubleVector w(this->distribution_pt(), 0.0);
  
     // Setup preconditioner only if we're not re-solving
     if (!Resolving)
     {
       // only setup the preconditioner before solve if require
       if (Setup_preconditioner_before_solve)
       {
         // Setup preconditioner from the Jacobian matrix
         double t_start_prec = TimingHelpers::timer();
  
         // do not setup
         preconditioner_pt()->setup(matrix_pt);
  
         // Doc time for setup of preconditioner
         double t_end_prec = TimingHelpers::timer();
         Preconditioner_setup_time = t_end_prec - t_start_prec;
  
         if (Doc_time)
         {
           oomph_info << "Time for setup of preconditioner  [sec]: "
                      << Preconditioner_setup_time << std::endl;
         }
       }
     }
     else
     {
       if (Doc_time)
       {
         oomph_info << "Setup of preconditioner is bypassed in resolve mode"
                    << std::endl;
       }
     }
  
     // solve b-Jx = Mr for r (assumes x = 0);
     DoubleVector r(this->distribution_pt(), 0.0);
     if (Preconditioner_LHS)
     {
       // Start the timer
       double t_start_prec = TimingHelpers::timer();
  
       // Apply the preconditioner
       preconditioner_pt()->preconditioner_solve(rhs, r);
  
       // Calculate the time taken for the preconditioner solve
       Preconditioner_application_time +=
         (TimingHelpers::timer() - t_start_prec);
     }
     else
     {
       r = rhs;
     }
     double normb = 0;
  
     // compute norm(r)
     double* r_pt = r.values_pt();
     for (unsigned i = 0; i < n_dof; i++)
     {
       normb += r_pt[i] * r_pt[i];
     }
     normb = sqrt(normb);
  
     // set beta (the initial residual)
     double beta = normb;
  
     // compute initial relative residual
     if (normb == 0.0) normb = 1;
     resid = beta / normb;
  
     // if required will document convergence history to screen or file (if
     // stream open)
     if (Doc_convergence_history)
     {
       if (!Output_file_stream.is_open())
       {
         oomph_info << 0 << " " << resid << std::endl;
       }
       else
       {
         Output_file_stream << 0 << " " << resid << std::endl;
       }
     }
  
     // if GMRES converges immediately
     if (resid <= Tolerance)
     {
       if (Doc_time)
       {
         oomph_info << "GMRES converged immediately. Normalised residual norm: "
                    << resid << std::endl;
       }
       // Doc time for solver
       double t_end = TimingHelpers::timer();
       Solution_time = t_end - t_start;
  
       if (Doc_time)
       {
         // Doc the time taken for the preconditioner applications
         oomph_info << "Time for all preconditioner applications [sec]: "
                    << Preconditioner_application_time
                    << "\n\nTime for solve with GMRES  [sec]: " << Solution_time
                    << std::endl;
       }
       return;
     }
  
  
     // initialise vector of orthogonal basis vectors (v) and upper hessenberg
     // matrix H
     // NOTE: for implementation purpose the upper hessenberg matrix indexes are
     // are swapped so the matrix is effectively transposed
     Vector<DoubleVector> v;
     v.resize(Restart + 1);
     Vector<Vector<double>> H(Restart + 1);
  
     // while...
     while (iter <= Max_iter)
     {
       // set zeroth basis vector v[0] to r/beta
       v[0].build(this->distribution_pt(), 0.0);
       double* v0_pt = v[0].values_pt();
       for (unsigned i = 0; i < n_dof; i++)
       {
         v0_pt[i] = r_pt[i] / beta;
       }
  
       //
       s[0] = beta;
  
       // inner iteration counter for restarted version
       unsigned iter_restart;
  
       // perform iterations
       for (iter_restart = 0; iter_restart < Restart && iter <= Max_iter;
            iter_restart++, iter++)
       {
         // resize next column of upper hessenberg matrix
         H[iter_restart].resize(iter_restart + 2);
  
         // solve Jv[i] = Mw for w
         {
           DoubleVector temp(this->distribution_pt(), 0.0);
           if (Preconditioner_LHS)
           {
             // solve Jv[i] = Mw for w
             matrix_pt->multiply(v[iter_restart], temp);
  
             // Start the timer
             double t_start_prec = TimingHelpers::timer();
  
             // Apply the preconditioner
             preconditioner_pt()->preconditioner_solve(temp, w);
  
             // Calculate the time taken for the preconditioner solve
             Preconditioner_application_time +=
               (TimingHelpers::timer() - t_start_prec);
           }
           else
           {
             // Start the timer
             double t_start_prec = TimingHelpers::timer();
  
             // w=JM^{-1}v by saad p270
             preconditioner_pt()->preconditioner_solve(v[iter_restart], temp);
  
             // Calculate the time taken for the preconditioner solve
             Preconditioner_application_time +=
               (TimingHelpers::timer() - t_start_prec);
  
             // Do a matrix multiplication
             matrix_pt->multiply(temp, w);
           }
         }
  
         //
         double* w_pt = w.values_pt();
         for (unsigned k = 0; k <= iter_restart; k++)
         {
           //
           H[iter_restart][k] = 0.0;
           double* vk_pt = v[k].values_pt();
           for (unsigned i = 0; i < n_dof; i++)
           {
             H[iter_restart][k] += w_pt[i] * vk_pt[i];
           }
  
           //
           for (unsigned i = 0; i < n_dof; i++)
           {
             w_pt[i] -= H[iter_restart][k] * vk_pt[i];
           }
         }
  
         //
         {
           double temp_norm_w = 0.0;
           for (unsigned i = 0; i < n_dof; i++)
           {
             temp_norm_w += w_pt[i] * w_pt[i];
           }
           temp_norm_w = sqrt(temp_norm_w);
           H[iter_restart][iter_restart + 1] = temp_norm_w;
         }
  
         //
         v[iter_restart + 1].build(this->distribution_pt(), 0.0);
         double* v_pt = v[iter_restart + 1].values_pt();
         for (unsigned i = 0; i < n_dof; i++)
         {
           v_pt[i] = w_pt[i] / H[iter_restart][iter_restart + 1];
         }
  
         //
         for (unsigned k = 0; k < iter_restart; k++)
         {
           apply_plane_rotation(
             H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
         }
         generate_plane_rotation(H[iter_restart][iter_restart],
                                 H[iter_restart][iter_restart + 1],
                                 cs[iter_restart],
                                 sn[iter_restart]);
         apply_plane_rotation(H[iter_restart][iter_restart],
                              H[iter_restart][iter_restart + 1],
                              cs[iter_restart],
                              sn[iter_restart]);
         apply_plane_rotation(s[iter_restart],
                              s[iter_restart + 1],
                              cs[iter_restart],
                              sn[iter_restart]);
  
         // compute current residual
         beta = std::fabs(s[iter_restart + 1]);
  
         // compute relative residual
         resid = beta / normb;
  
         // if required will document convergence history to screen or file (if
         // stream open)
         if (Doc_convergence_history)
         {
           if (!Output_file_stream.is_open())
           {
             oomph_info << iter << " " << resid << std::endl;
           }
           else
           {
             Output_file_stream << iter << " " << resid << std::endl;
           }
         }
  
         // if required tolerance found
         if (resid < Tolerance)
         {
           // update result vector
           update(iter_restart, H, s, v, solution);
  
           // document convergence
           if (Doc_time)
           {
             oomph_info << std::endl;
             oomph_info << "GMRES converged (1). Normalised residual norm: "
                        << resid << std::endl;
             oomph_info << "Number of iterations to convergence: " << iter
                        << std::endl;
             oomph_info << std::endl;
           }
  
           // Doc time for solver
           double t_end = TimingHelpers::timer();
           Solution_time = t_end - t_start;
  
           Iterations = iter;
  
           if (Doc_time)
           {
             // Doc the time taken for the preconditioner applications
             oomph_info << "Time for all preconditioner applications [sec]: "
                        << Preconditioner_application_time
                        << "\n\nTime for solve with GMRES  [sec]: "
                        << Solution_time << std::endl;
           }
           return;
         }
       }
  
       // update
       if (iter_restart > 0) update((iter_restart - 1), H, s, v, solution);
  
       // solve Mr = (b-Jx) for r
       {
         DoubleVector temp(this->distribution_pt(), 0.0);
         matrix_pt->multiply(solution, temp);
         double* temp_pt = temp.values_pt();
         const double* rhs_pt = rhs.values_pt();
         for (unsigned i = 0; i < n_dof; i++)
         {
           temp_pt[i] = rhs_pt[i] - temp_pt[i];
         }
  
         if (Preconditioner_LHS)
         {
           // Start the timer
           double t_start_prec = TimingHelpers::timer();
  
           preconditioner_pt()->preconditioner_solve(temp, r);
  
           // Calculate the time taken for the preconditioner solve
           Preconditioner_application_time +=
             (TimingHelpers::timer() - t_start_prec);
         }
       }
  
       // compute current residual
       beta = 0.0;
       r_pt = r.values_pt();
       for (unsigned i = 0; i < n_dof; i++)
       {
         beta += r_pt[i] * r_pt[i];
       }
       beta = sqrt(beta);
  
       // if relative residual within tolerance
       resid = beta / normb;
       if (resid < Tolerance)
       {
         if (Doc_time)
         {
           oomph_info << std::endl;
           oomph_info << "GMRES converged (2). Normalised residual norm: "
                      << resid << std::endl;
           oomph_info << "Number of iterations to convergence: " << iter
                      << std::endl;
           oomph_info << std::endl;
         }
  
         // Doc time for solver
         double t_end = TimingHelpers::timer();
         Solution_time = t_end - t_start;
  
         Iterations = iter;
  
         if (Doc_time)
         {
           // Doc the time taken for the preconditioner applications
           oomph_info << "Time for all preconditioner applications [sec]: "
                      << Preconditioner_application_time
                      << "\n\nTime for solve with GMRES  [sec]: "
                      << Solution_time << std::endl;
         }
         return;
       }
     }
  
     // DRAIG: I want to delete this but GMRES doesn't compute the residuals
     // properly after this scope when it doesn't converge within the specified
     // number of iterations. Not sure why. If you're bored, try to figure
     // out why.
     {
       DoubleVector temp(this->distribution_pt(), 0.0);
       matrix_pt->multiply(solution, temp);
       temp *= -1.0;
       temp += rhs;
       resid = temp.norm() / normb;
     }
  
     // otherwise GMRES failed convergence
     oomph_info << std::endl;
     oomph_info << "GMRES did not converge to required tolerance! " << std::endl;
     oomph_info << "Returning with normalised residual norm: " << resid
                << std::endl;
     oomph_info << "after " << Max_iter << " iterations." << std::endl;
     oomph_info << std::endl;
  
  
     if (Throw_error_after_max_iter)
     {
       std::string err = "Solver failed to converge and you requested an error";
       err += " on convergence failures.";
       throw OomphLibError(
         err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
     }
  
     return;
  
   } // End GMRES

References beta, oomph::DoubleVector::built(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), boost::multiprecision::fabs(), H, i, Global_Variables::Iterations, k, oomph::BlackBoxFDNewtonSolver::Max_iter, oomph::DoubleMatrixBase::multiply(), oomph::DoubleVector::norm(), oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, UniformPSDSelfTest::r, s, BiharmonicTestFunctions1::solution(), sqrt(), oomph::Global_string_for_annotation::string(), oomph::TimingHelpers::timer(), v, oomph::DoubleVector::values_pt(), and w.

Referenced by oomph::GMRES< MATRIX >::solve().

◆ update()

template<typename MATRIX >

void oomph::GMRES< MATRIX >::update	(	const unsigned &	k,
		const Vector< Vector< double >> &	H,
		const Vector< double > &	s,
		const Vector< DoubleVector > &	v,
		DoubleVector &	x
	)

inlineprivate

Helper function to update the result vector using the result, x=x_0+V_m*y

     {
       // Make a local copy of s
       Vector<double> y(s);
  
       // Backsolve:
       for (int i = int(k); i >= 0; i--)
       {
         // Divide the i-th entry of y by the i-th diagonal entry of H
         y[i] /= H[i][i];
  
         // Loop over the previous values of y and update
         for (int j = i - 1; j >= 0; j--)
         {
           // Update the j-th entry of y
           y[j] -= H[i][j] * y[i];
         }
       } // for (int i=int(k);i>=0;i--)
  
       // Store the number of rows in the result vector
       unsigned n_x = x.nrow();
  
       // Build a temporary vector with entries initialised to 0.0
       DoubleVector temp(x.distribution_pt(), 0.0);
  
       // Build a temporary vector with entries initialised to 0.0
       DoubleVector z(x.distribution_pt(), 0.0);
  
       // Get access to the underlying values
       double* temp_pt = temp.values_pt();
  
       // Calculate x=Vy
       for (unsigned j = 0; j <= k; j++)
       {
         // Get access to j-th column of v
         const double* vj_pt = v[j].values_pt();
  
         // Loop over the entries of the vector, temp
         for (unsigned i = 0; i < n_x; i++)
         {
           temp_pt[i] += vj_pt[i] * y[j];
         }
       } // for (unsigned j=0;j<=k;j++)
  
       // If we're using LHS preconditioning
       if (Preconditioner_LHS)
       {
         // Since we're using LHS preconditioning the preconditioner is applied
         // to the matrix and RHS vector so we simply update the value of x
         x += temp;
       }
       // If we're using RHS preconditioning
       else
       {
         // Start the timer
         double t_start_prec = TimingHelpers::timer();
  
         // Since we're using RHS preconditioning the preconditioner is applied
         // to the solution vector
         preconditioner_pt()->preconditioner_solve(temp, z);
  
         // Calculate the time taken for the preconditioner solve
         Preconditioner_application_time +=
           (TimingHelpers::timer() - t_start_prec);
  
         // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
         // sparse linear systems", p.284]
         x += z;
       }
     } // End of update

private

The number of iterations before the iteration proceedure is restarted if iteration restart is used

The documentation for this class was generated from the following files:

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

template<typename MATRIX> class oomph::GMRES< MATRIX >

Constructor & Destructor Documentation

◆ GMRES() [1/2]

◆ ~GMRES()

◆ GMRES() [2/2]

Member Function Documentation

◆ apply_plane_rotation()

◆ clean_up_memory()

◆ disable_iteration_restart()

◆ disable_resolve()

◆ enable_computation_of_gradient()

◆ enable_iteration_restart()

◆ generate_plane_rotation()

◆ iteration_restart()

◆ iterations()

◆ operator=()

◆ resolve()

◆ set_preconditioner_LHS()

◆ set_preconditioner_RHS()

◆ solve() [1/3]

◆ solve() [2/3]

◆ solve() [3/3]

◆ solve_helper()

◆ update()

Member Data Documentation

◆ Iteration_restart

◆ Iterations

◆ Matrix_can_be_deleted

◆ Matrix_pt

◆ Preconditioner_application_time

◆ Preconditioner_LHS

◆ Resolving

◆ Restart

template<typename MATRIX>
class oomph::GMRES< MATRIX >