The GMRES method. More...

#include <iterative_linear_solver.h>

Inheritance diagram for oomph::AugmentedProblemGMRES:

Public Member Functions
	AugmentedProblemGMRES (DoubleVector b_pt, DoubleVector c_pt, double x_pt, double rhs_pt)
	Constructor. More...

virtual	~AugmentedProblemGMRES ()
	Destructor: Clean up storage. More...

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

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

void	disable_resolve ()
	Overload disable resolve so that it cleans up memory too. 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)

virtual void	enable_computation_of_gradient ()

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	clean_up_memory ()
	Cleanup data that's stored for resolve (if any has been stored) More...

void	augmented_matrix_multiply (CRDoubleMatrix *matrix_pt, const DoubleVector &x, DoubleVector &soln)
	Multiply the vector x by the augmented system matrix. More...

void	apply_schur_complement_preconditioner (const DoubleVector &rhs, DoubleVector &soln)

void	solve_helper (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
	General interface to solve function. 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...

CRDoubleMatrix *	Matrix_pt
	Pointer to matrix. More...

DoubleVector *	B_pt
	Pointer to the column vector in the bordered system. More...

DoubleVector *	C_pt
	Pointer to the row vector in the bordered system. More...

double *	X_pt
	Pointer to the last entry of the LHS vector in the bordered system. More...

double *	Rhs_pt
	Pointer to the last entry of the RHS vector in the bordered system. More...

double	Schur_complement_scalar
	The scalar component of the Schur complement preconditioner. More...

bool	Resolving

bool	Matrix_can_be_deleted

bool	Preconditioner_LHS

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

The GMRES method.

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

Constructor & Destructor Documentation

◆ AugmentedProblemGMRES() [1/2]

oomph::AugmentedProblemGMRES::AugmentedProblemGMRES	(	DoubleVector *	b_pt,
		DoubleVector *	c_pt,
		double *	x_pt,
		double *	rhs_pt
	)

inline

Constructor.

       : Iterations(0),
         Matrix_pt(0),
         B_pt(b_pt),
         C_pt(c_pt),
         X_pt(x_pt),
         Rhs_pt(rhs_pt),
         Schur_complement_scalar(1.0),
         Resolving(false),
         Matrix_can_be_deleted(true)
     {
       Preconditioner_LHS = true;
       Iteration_restart = false;
     }

References Iteration_restart, and Preconditioner_LHS.

◆ ~AugmentedProblemGMRES()

virtual oomph::AugmentedProblemGMRES::~AugmentedProblemGMRES ( )

inlinevirtual

Destructor: Clean up storage.

     {
       clean_up_memory();
     }

References clean_up_memory().

◆ AugmentedProblemGMRES() [2/2]

oomph::AugmentedProblemGMRES::AugmentedProblemGMRES ( const AugmentedProblemGMRES & )

delete

Broken copy constructor.

Member Function Documentation

◆ apply_plane_rotation()

void oomph::AugmentedProblemGMRES::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;
     }

Referenced by solve_helper().

◆ apply_schur_complement_preconditioner()

void oomph::AugmentedProblemGMRES::apply_schur_complement_preconditioner	(	const DoubleVector &	rhs,
		DoubleVector &	soln
	)

inlineprivate

Apply the block-diagonal Schur complement preconditioner to compute the LHS which has size N+1 (the augmented system size)

     {
       // How many dofs are there in the non-augmented system?
       unsigned n_dof = this->distribution_pt()->nrow();
  
       // Compute the final entry of r first
       soln[n_dof] = rhs[n_dof] / Schur_complement_scalar;
  
       // Do we want to use a block-diagonal approximation? The default is to
       // use a block upper-triangular approximation for the preconditioner
       bool use_block_diagonal_preconditioner = false;
  
       // Allocate space for (the non-augmented part of) r satisfying b-Jx = Mr
       DoubleVector r_small(this->distribution_pt(), 0.0);
  
       // Allocate space for the non-augmented system part of the RHS vector
       DoubleVector rhs_small(this->distribution_pt(), 0.0);
  
       // Loop over the first n_dof entries of the RHS vector
       for (unsigned i = 0; i < n_dof; i++)
       {
         // If we want to use the block-diagonal preconditioner
         if (use_block_diagonal_preconditioner)
         {
           // Copy the i-th entry over
           rhs_small[i] = rhs[i];
         }
         // If we want to use the block upper-triangular preconditioner
         else
         {
           // Copy the i-th entry over
           rhs_small[i] = rhs[i] - soln[n_dof] * (*B_pt)[i];
         }
       } // for (unsigned i=0;i<n_dof;i++)
  
       // Apply the preconditioner
       preconditioner_pt()->preconditioner_solve(rhs_small, r_small);
  
       // Loop over the first n_dof entries of the solution vector
       for (unsigned i = 0; i < n_dof; i++)
       {
         // Copy the i-th entry over
         soln[i] = r_small[i];
       }
     } // End of apply_schur_complement_preconditioner

References oomph::DistributableLinearAlgebraObject::distribution_pt(), i, oomph::LinearAlgebraDistribution::nrow(), oomph::IterativeLinearSolver::preconditioner_pt(), oomph::Preconditioner::preconditioner_solve(), and Schur_complement_scalar.

Referenced by solve_helper(), and update().

◆ augmented_matrix_multiply()

void oomph::AugmentedProblemGMRES::augmented_matrix_multiply	(	CRDoubleMatrix *	matrix_pt,
		const DoubleVector &	x,
		DoubleVector &	soln
	)

inlineprivate

Multiply the vector x by the augmented system matrix.

     {
       // How many dofs are there in the non-augmented system?
       unsigned n_dof = matrix_pt->nrow();
  
       // Allocate space for the non-augmented component of x
       DoubleVector x_small(this->distribution_pt(), 0.0);
  
       // Loop over the first n_dof entries of x
       for (unsigned i = 0; i < n_dof; i++)
       {
         // Copy the i-th entry over
         x_small[i] = x[i];
       }
  
       // Allocate space for the product of the matrix and x_small
       DoubleVector a_prod_xsmall(this->distribution_pt(), 0.0);
  
       // Now multiply the matrix and x_small
       matrix_pt->multiply(x_small, a_prod_xsmall);
  
       // Get the scalar component of x associated with the system augmentation
       double y = x[n_dof];
  
       // Loop over the first n_dof entries of soln
       for (unsigned i = 0; i < n_dof; i++)
       {
         // Compute the i-th entry
         soln[i] = a_prod_xsmall[i] + y * (*B_pt)[i];
       }
  
       // Compute the final entry of soln
       soln[n_dof] = C_pt->dot(x_small);
     } // End of augmented_matrix_multiply

References C_pt, oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::DoubleVector::dot(), i, oomph::CRDoubleMatrix::multiply(), oomph::CRDoubleMatrix::nrow(), plotDoE::x, and y.

Referenced by solve_helper().

◆ clean_up_memory()

void oomph::AugmentedProblemGMRES::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;
       }
     } // End of clean_up_memory

References Matrix_can_be_deleted, and Matrix_pt.

Referenced by disable_resolve(), solve(), and ~AugmentedProblemGMRES().

◆ disable_iteration_restart()

void oomph::AugmentedProblemGMRES::disable_iteration_restart ( )

inline

Switches off iteration restart.

     {
       Iteration_restart = false;
     } // End of disable_iteration_restart

References Iteration_restart.

◆ disable_resolve()

void oomph::AugmentedProblemGMRES::disable_resolve ( )

inlinevirtual

Overload disable resolve so that it cleans up memory too.

Reimplemented from oomph::LinearSolver.

     {
       LinearSolver::disable_resolve();
       clean_up_memory();
     } // End of disable_resolve

References clean_up_memory(), and oomph::LinearSolver::disable_resolve().

◆ enable_iteration_restart()

void oomph::AugmentedProblemGMRES::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 Iteration_restart.

◆ generate_plane_rotation()

void oomph::AugmentedProblemGMRES::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().

Referenced by solve_helper().

◆ iteration_restart()

bool oomph::AugmentedProblemGMRES::iteration_restart ( ) const

inline

access function indicating whether restarted GMRES is used

     {
       return Iteration_restart;
     }

References Iteration_restart.

◆ iterations()

unsigned oomph::AugmentedProblemGMRES::iterations ( ) const

inlinevirtual

Number of iterations taken.

Implements oomph::IterativeLinearSolver.

     {
       return Iterations;
     }

References Iterations.

◆ operator=()

void oomph::AugmentedProblemGMRES::operator= ( const AugmentedProblemGMRES & )

delete

Broken assignment operator.

◆ resolve()

void oomph::AugmentedProblemGMRES::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 Matrix_pt, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Resolving, and solve().

◆ set_preconditioner_LHS()

void oomph::AugmentedProblemGMRES::set_preconditioner_LHS ( )

inline

Set left preconditioning (the default)

     {
       Preconditioner_LHS = true;
     }

References Preconditioner_LHS.

◆ set_preconditioner_RHS()

void oomph::AugmentedProblemGMRES::set_preconditioner_RHS ( )

inline

Enable right preconditioning.

     {
       Preconditioner_LHS = false;
     }

References Preconditioner_LHS.

◆ solve() [1/3]

void oomph::AugmentedProblemGMRES::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.

     {
       // Upcast to a CRDoubleMatrix
       CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
  
       // If we can't upcast to a CRDoubleMatrix then this won't work...
       if (cr_matrix_pt == 0)
       {
         // Throw an error
         throw OomphLibError("Can't upcast input matrix to a CRDoubleMatrix!",
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
  
       // Set up the distribution
       this->build_distribution(cr_matrix_pt->distribution_pt());
  
       // Store the matrix if required
       if ((Enable_resolve) && (!Resolving))
       {
         // Store a pointer to the upcasted matrix
         Matrix_pt = cr_matrix_pt;
  
         // Matrix has been passed in from the outside so we must not
         // delete it
         Matrix_can_be_deleted = false;
       }
  
       // Call the helper function
       this->solve_helper(matrix_pt, rhs, solution);
     } // End of solve

References oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Enable_resolve, Matrix_can_be_deleted, Matrix_pt, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Resolving, BiharmonicTestFunctions1::solution(), and solve_helper().

◆ solve() [2/3]

void oomph::AugmentedProblemGMRES::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]

void oomph::AugmentedProblemGMRES::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.

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// 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 (in the non-augmented problem)
     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();
  
     // Create a distribution
     LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
  
     // ...now build it
     this->build_distribution(dist);
  
     // Get the Jacobian matrix and the nonlinear residual vector
     Matrix_pt = new CRDoubleMatrix;
     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;
     }
  
     // Reset the Schur complement scalar entry value
     Schur_complement_scalar = 1.0;
  
     // Call linear algebra-style solver. If the result distribution is wrong,
     // then redistribute before the solve and return to original distribution
     // afterwards
     if ((!(*result.distribution_pt() == *this->distribution_pt())) &&
         result.built())
     {
       // Create a temporary copy of the current result distribution
       LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
  
       // Re-build the result vector as an augmented vector
       result.build(dist, 0.0);
  
       // Call the auxilliary helper function to do the actual solve
       this->solve_helper(Matrix_pt, f, result);
  
       // Re-distribute result vector
       result.redistribute(&temp_global_dist);
     }
     // Otherwise just solve
     else
     {
       // Call the auxilliary helper function to do the actual solve
       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::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::built(), clean_up_memory(), oomph::Problem::communicator_pt(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Doc_time, oomph::LinearSolver::Enable_resolve, f(), oomph::Problem::get_jacobian(), oomph::IterativeLinearSolver::Jacobian_setup_time, Matrix_can_be_deleted, Matrix_pt, oomph::Problem::ndof(), oomph::oomph_info, oomph::DoubleVector::redistribute(), Resolving, Schur_complement_scalar, solve_helper(), and oomph::TimingHelpers::timer().

Referenced by resolve().

◆ solve_helper()

void oomph::AugmentedProblemGMRES::solve_helper	(	DoubleMatrixBase *const &	input_matrix_pt,
		const DoubleVector &	rhs,
		DoubleVector &	lhs
	)

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 (in the non-augmented problem)
     unsigned n_dof = input_matrix_pt->nrow();
  
     // How many dofs in the augmented problem?
     unsigned n_aug_dof = n_dof + 1;
  
     // Get a handle to this Problems communicator
     OomphCommunicator* comm_pt = this->distribution_pt()->communicator_pt();
  
     // Set up the (non-augmented) distribution
     LinearAlgebraDistribution dist(comm_pt, n_dof, false);
  
     // Set up the (augmented) distribution
     LinearAlgebraDistribution aug_dist(comm_pt, n_aug_dof, false);
  
     // If the input vector doesn't have size N+1 then something has gone wrong
     if ((rhs.nrow() != n_dof) || (lhs.nrow() != n_dof))
     {
       // Create an output stream
       std::ostringstream error_message_stream;
  
       // Create the error message
       error_message_stream << "RHS vector has " << rhs.nrow()
                            << " rows and the "
                            << "LHS vector has " << rhs.nrow()
                            << " rows but\nthey "
                            << "should both have " << n_dof << " rows!"
                            << std::endl;
  
       // Throw an error
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
 #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*>(input_matrix_pt) != 0)
     {
       if (dynamic_cast<DistributableLinearAlgebraObject*>(input_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 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 (lhs.built())
     {
       if (!(*rhs.distribution_pt() == *lhs.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 (!lhs.built())
     {
       // Use the augmented distribution which has size N+1 (not the solvers
       // distribution which has size N)
       lhs.build(dist, 0.0);
     }
     // Otherwise initialise to zero
     else
     {
       lhs.initialise(0.0);
     }
  
     // Create a vector to store the augmented LHS vector entries
     DoubleVector solution(aug_dist, 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_aug_dof
     if (!Iteration_restart)
     {
       Restart = n_aug_dof;
     }
  
     // initialise vectors
     Vector<double> s(Restart + 1, 0);
     Vector<double> cs(Restart + 1);
     Vector<double> sn(Restart + 1);
     DoubleVector w(aug_dist, 0.0);
  
     // Set up the 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(input_matrix_pt);
  
         // Doc time for setup of preconditioner
         double t_end_prec = TimingHelpers::timer();
         Preconditioner_setup_time = t_end_prec - t_start_prec;
  
         // If we're meant to document timings
         if (Doc_time)
         {
           oomph_info << "Time for setup of preconditioner [sec]: "
                      << Preconditioner_setup_time << std::endl;
         }
       } // if (Setup_preconditioner_before_solve)
     }
     else
     {
       // If we're meant to document timings
       if (Doc_time)
       {
         oomph_info << "Setup of preconditioner is bypassed in resolve mode"
                    << std::endl;
       }
     } // if (!Resolving)
  
     // The default system matrix is simply the one that was passed in
     CRDoubleMatrix* matrix_pt = dynamic_cast<CRDoubleMatrix*>(input_matrix_pt);
  
     // Storage for the time taken for all preconditioner applications
     double t_prec_application_total = 0.0;
  
     //---------------------------------------------------------------
     // Compute the scalar component of the diagonal Schur complement
     // preconditioner. NOTE: Only compute it if it is currently set
     // to the default value of unity. It's quite unlikely that it
     // will actually have value one (up to machine precision) in
     // practice. This check ensures that we only compute the value
     // during the solve function and not during the resolve.
     //---------------------------------------------------------------
     // If the Schur complement value has its default value of one
     if (std::fabs(Schur_complement_scalar - 1.0) < 1.0e-14)
     {
       // If the bordering vectors haven't actually been properly created
       if ((B_pt == 0) || (C_pt == 0))
       {
         // Throw an error
         throw OomphLibError("Bordering vectors have not been created!",
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
       // Make sure the bordering vectors have the right size
       else if ((B_pt->nrow() != n_dof) || (C_pt->nrow() != n_dof))
       {
         // Throw an error
         throw OomphLibError("Bordering vectors do not have the right size!",
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
  
       // Allocate storage for A^{-1}b
       DoubleVector a_inv_b(this->distribution_pt(), 0.0);
  
       // Compute A^{-1}b with the preconditioner approximation to A^{-1}
       preconditioner_pt()->preconditioner_solve(*B_pt, a_inv_b);
  
       // Now compute the scalar component of the Schur complement preconditioner
       Schur_complement_scalar = C_pt->dot(a_inv_b);
  
       // If the Schur complement scalar value is pratically zero or is NaN
       if ((std::fabs(Schur_complement_scalar) < 1.0e-14) ||
           std::isnan(Schur_complement_scalar))
       {
         // Reset it to have value 2.0 (not 1.0 otherwise we would end up
         // recomputing it every resolve) as we "should" be able to use
         // any scalar value and end up with the same number of distinct
         // eigenvalues for the preconditioned system).
         Schur_complement_scalar = 2.0;
       }
     } // if (std::fabs(Schur_complement_scalar-1.0)<1.0e-14)
  
     // Storage for the augmented RHS vector
     DoubleVector rhs_aug(aug_dist, 0.0);
  
     // Loop over the first n_dof entries of the RHS vector
     for (unsigned i = 0; i < n_dof; i++)
     {
       // Copy the i-th entry over
       rhs_aug[i] = rhs[i];
     }
  
     // Update the RHS entry associated with the augmented equation
     rhs_aug[n_dof] = (*Rhs_pt);
  
     // Storage for b-Jx=Mr for r (assumes x = 0);
     DoubleVector r(aug_dist, 0.0);
  
     // If we're using LHS preconditioning
     if (Preconditioner_LHS)
     {
       // Initialise timer
       double t_prec_application_start = TimingHelpers::timer();
  
       // Use the helper function to apply the Schur complement preconditioner
       apply_schur_complement_preconditioner(rhs_aug, r);
  
       // Doc time for setup of preconditioner
       double t_prec_application_end = TimingHelpers::timer();
  
       // Update the preconditioner application time total
       t_prec_application_total +=
         (t_prec_application_end - t_prec_application_start);
     }
     // If we're using RHS preconditioning
     else
     {
       // We don't need to do any preconditioning here; just copy the vector over
       r = rhs_aug;
     } // if (Preconditioner_LHS)
  
     // Compute the norm of the RHS vector
     double normb = r.norm();
  
     // Set beta (the initial residual)
     double beta = normb;
  
     // If the RHS vector is (numerically equivalent to) the zero vector
     if (normb == 0.0)
     {
       // Initialise the norm of b as unity
       normb = 1.0;
     }
  
     // Compute the 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 << 0 << " " << resid << std::endl;
       }
       else
       {
         Output_file_stream << 0 << " " << resid << std::endl;
       }
     } // if (Doc_convergence_history)
  
     // If GMRES converges immediately
     if (resid <= Tolerance)
     {
       if (Doc_time)
       {
         oomph_info << "AugmentedProblemGMRES 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 all preconditioner applications [sec]: "
           << t_prec_application_total
           << "\nTotal time for solve with AugmentedProblemGMRES [sec]: "
           << Solution_time << std::endl;
       }
  
       // We're done so finish here
       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)
     {
       // Copy the values of r into the zeroth basis vector v[0]
       v[0] = r;
  
       // Now scale it so the zeroth basis vector v[0] is r/beta
       v[0] /= beta;
  
       // Set the first entry of s
       s[0] = beta;
  
       // Counter for the inner iteration counter (restarted version)
       unsigned iter_restart = 0;
  
       // Start performing 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
         {
           // Allocate space for the product J*v[i]
           DoubleVector temp(aug_dist, 0.0);
  
           // If we're using LHS preconditioning
           if (Preconditioner_LHS)
           {
             // Compute J*v[i] where J is now the augmented system matrix
             augmented_matrix_multiply(matrix_pt, v[iter_restart], temp);
  
             // Initialise timer
             double t_prec_application_start = TimingHelpers::timer();
  
             // Use the helper function to apply the Schur complement
             // preconditioner
             apply_schur_complement_preconditioner(temp, w);
  
             // End timer
             double t_prec_application_end = TimingHelpers::timer();
  
             // Update the preconditioner application time total
             t_prec_application_total +=
               (t_prec_application_end - t_prec_application_start);
           }
           else
           {
             // Initialise timer
             double t_prec_application_start = TimingHelpers::timer();
  
             // Compute w=JM^{-1}v by saad p270
             apply_schur_complement_preconditioner(v[iter_restart], temp);
  
             // Doc time for setup of preconditioner
             double t_prec_application_end = TimingHelpers::timer();
  
             // Update the preconditioner application time total
             t_prec_application_total +=
               (t_prec_application_end - t_prec_application_start);
  
             // Calculate w=Jv_m where v_m=M^{-1}v
             augmented_matrix_multiply(matrix_pt, temp, w);
           }
         } // Solve Jv[i]=Mw for w
  
         // Get a pointer to the entries of w
         double* w_pt = w.values_pt();
  
         // Loop over the rows of the Hessenberg matrix
         for (unsigned k = 0; k <= iter_restart; k++)
         {
           // Compute the (k,iter_restart)-th entry of the Hessenberg matrix
           // (which is stored in its transposed form)
           H[iter_restart][k] = w.dot(v[k]);
  
           // Get a pointer to the entries of the k-th basis vector
           double* vk_pt = v[k].values_pt();
  
           // Loop over the entries of w
           for (unsigned i = 0; i < n_aug_dof; i++)
           {
             // Update the i-th entry of w
             w_pt[i] -= H[iter_restart][k] * vk_pt[i];
           }
         } // for (unsigned k=0;k<=iter_restart;k++)
  
         // Calculate the subdiagonal Hessenberg entry
         H[iter_restart][iter_restart + 1] = w.norm();
  
         // Copy the entries of w into the next basis vector
         v[iter_restart + 1] = w;
  
         // Now normalise the basis vector with the appropriate Hessenberg entry
         v[iter_restart + 1] /= H[iter_restart][iter_restart + 1];
  
         // Loop over the rows of the Hessenberg matrix
         for (unsigned k = 0; k < iter_restart; k++)
         {
           // Apply a Givens rotation to the appropriate entries of H
           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 the current residual
         beta = std::fabs(s[iter_restart + 1]);
  
         // compute relative residual
         resid = beta / normb;
  
         // DRAIG: Delete...
         // oomph_info << "Residual at iteration "
         //       << iter << ": " << resid << std::endl;
  
         // 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 (Doc_convergence_history)
  
         // If required tolerance found
         if (resid < Tolerance)
         {
           // update result vector
           update(iter_restart, H, s, v, solution);
  
           // Loop over the first n_dof entries of the RHS vector
           for (unsigned i = 0; i < n_dof; i++)
           {
             // Copy the i-th entry over
             lhs[i] = solution[i];
           }
  
           // Update the bifurcation parameter
           (*X_pt) -= solution[n_dof];
  
           // If we're meant to document convergence
           if (Doc_time)
           {
             // Output the convergence to screen
             oomph_info << "\nAugmentedProblemGMRES converged (1). Normalised "
                        << "residual norm: " << resid
                        << "\nNumber of iterations "
                        << "to convergence: " << iter << "\n"
                        << std::endl;
           }
  
           // Doc time for solver
           double t_end = TimingHelpers::timer();
           Solution_time = t_end - t_start;
  
           Iterations = iter;
  
           // If we're meant to document convergence
           if (Doc_time)
           {
             oomph_info
               << "Time for all preconditioner applications [sec]: "
               << t_prec_application_total
               << "\nTotal time for solve with AugmentedProblemGMRES  [sec]: "
               << Solution_time << std::endl;
           }
  
           // We're done now so jump out here
           return;
         }
 #ifdef PARANOID
         // If the tolerance wasn't met
         else
         {
           // Warn the user if the sub-diagonal Hessenberg entry was almost
           // zero...
           if (std::fabs(H[iter_restart][iter_restart + 1]) < 1.0e-12)
           {
             // Create an output stream
             std::ostringstream warning_message_stream;
  
             // Create the warning message
             warning_message_stream << "Hessenberg subdiagonal entry (="
                                    << H[iter_restart][iter_restart + 1]
                                    << ") almost zero at iteration "
                                    << iter_restart << "\nbut residual "
                                    << "is only " << resid << " which does "
                                    << "not meet the tolerance (" << Tolerance
                                    << ")!" << std::endl;
  
             // Warn the user
             OomphLibWarning(warning_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
           }
         } // if (resid<Tolerance)
 #endif
       } // for (iter_restart=0;iter_restart<Restart&&iter<=Max_iter;...
  
       // If we're not on the first iteration
       if (iter_restart > 0)
       {
         // Do an update
         update((iter_restart - 1), H, s, v, solution);
       }
  
       // Solve Mr = (b-Jx) for r
       {
         DoubleVector temp(aug_dist, 0.0);
         augmented_matrix_multiply(matrix_pt, solution, temp);
         double* temp_pt = temp.values_pt();
         const double* rhs_pt = rhs.values_pt();
         for (unsigned i = 0; i < n_dof; i++)
         {
           // Update the i-th entry of temp to be b-Jx
           temp_pt[i] = rhs_pt[i] - temp_pt[i];
         }
  
         // If we're using LHS preconditioning
         if (Preconditioner_LHS)
         {
           // Initialise timer
           double t_prec_application_start = TimingHelpers::timer();
  
           // Use the helper function to apply the Schur complement
           // preconditioner
           apply_schur_complement_preconditioner(temp, r);
  
           // Doc time for setup of preconditioner
           double t_prec_application_end = TimingHelpers::timer();
  
           // Update the preconditioner application time total
           t_prec_application_total +=
             (t_prec_application_end - t_prec_application_start);
         }
       } // Solve Mr = (b-Jx) for r
  
       // Compute the current residual
       beta = r.norm();
  
       // if relative residual within tolerance
       resid = beta / normb;
       if (resid < Tolerance)
       {
         // If we're meant to document convergence
         if (Doc_time)
         {
           oomph_info << "\nAugmentedProblemGMRES converged (2). Normalised "
                      << "residual norm: " << resid << "\nNumber of iterations "
                      << "to convergence: " << iter << "\n"
                      << std::endl;
         }
  
         // Loop over the first n_dof entries of the RHS vector
         for (unsigned i = 0; i < n_dof; i++)
         {
           // Copy the i-th entry over
           lhs[i] = solution[i];
         }
  
         // Update the bifurcation parameter
         (*X_pt) -= solution[n_dof];
  
         // Doc time for solver
         double t_end = TimingHelpers::timer();
         Solution_time = t_end - t_start;
  
         Iterations = iter;
  
         // If we're meant to document convergence
         if (Doc_time)
         {
           oomph_info
             << "Time for all preconditioner applications [sec]: "
             << t_prec_application_total
             << "\nTotal time for solve with AugmentedProblemGMRES [sec]: "
             << Solution_time << std::endl;
         }
  
         return;
       } // if (resid<Tolerance)
     } // while (iter<=Max_iter)
  
     // DRAIG: Delete...
     {
       DoubleVector temp(aug_dist, 0.0);
       augmented_matrix_multiply(matrix_pt, solution, temp);
       temp *= -1.0;
       temp += rhs;
       resid = temp.norm() / normb;
     }
  
     // Loop over the first n_dof entries of the RHS vector
     for (unsigned i = 0; i < n_dof; i++)
     {
       // Copy the i-th entry over
       lhs[i] = solution[i];
     }
  
     // Update the bifurcation parameter
     (*X_pt) -= solution[n_dof];
  
     // Otherwise AugmentedProblemGMRES failed convergence
     oomph_info << "\nAugmentedProblemGMRES did not converge to required "
                << "tolerance!\nReturning with normalised residual norm: "
                << resid << "\nafter " << Max_iter << " iterations.\n"
                << std::endl;
  
     if (Throw_error_after_max_iter)
     {
       std::string error_message_string = "Solver failed to converge and you ";
       error_message_string += "requested an error on convergence failures.";
       throw OomphLibError(
         error_message_string, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
     }
  
     return;
   } // End AugmentedProblemGMRES

References apply_plane_rotation(), apply_schur_complement_preconditioner(), augmented_matrix_multiply(), B_pt, beta, oomph::DoubleVector::build(), oomph::DoubleVector::built(), C_pt, oomph::LinearAlgebraDistribution::communicator_pt(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::IterativeLinearSolver::Doc_convergence_history, oomph::LinearSolver::Doc_time, oomph::DoubleVector::dot(), boost::multiprecision::fabs(), generate_plane_rotation(), H, i, oomph::DoubleVector::initialise(), isnan, Iteration_restart, Iterations, k, oomph::IterativeLinearSolver::Max_iter, oomph::DoubleVector::norm(), oomph::DistributableLinearAlgebraObject::nrow(), oomph::DoubleMatrixBase::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::IterativeLinearSolver::Output_file_stream, Preconditioner_LHS, oomph::IterativeLinearSolver::preconditioner_pt(), oomph::IterativeLinearSolver::Preconditioner_setup_time, oomph::Preconditioner::preconditioner_solve(), UniformPSDSelfTest::r, Resolving, s, Schur_complement_scalar, oomph::Preconditioner::setup(), oomph::IterativeLinearSolver::Setup_preconditioner_before_solve, BiharmonicTestFunctions1::solution(), oomph::IterativeLinearSolver::Solution_time, oomph::Global_string_for_annotation::string(), oomph::IterativeLinearSolver::Throw_error_after_max_iter, oomph::TimingHelpers::timer(), oomph::IterativeLinearSolver::Tolerance, update(), v, oomph::DoubleVector::values_pt(), and w.

Referenced by solve().

◆ update()

void oomph::AugmentedProblemGMRES::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
       {
         // Since we're using RHS preconditioning the preconditioner is applied
         // to the solution vector
         apply_schur_complement_preconditioner(temp, z);
  
         // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
         // sparse linear systems", p.284]
         x += z;
       }
     } // End of update

References apply_schur_complement_preconditioner(), H, i, j, k, Preconditioner_LHS, s, v, oomph::DoubleVector::values_pt(), plotDoE::x, and y.

Referenced by smc.smc::recursiveBayesian(), and solve_helper().

Member Data Documentation

◆ B_pt

DoubleVector* oomph::AugmentedProblemGMRES::B_pt

private

Pointer to the column vector in the bordered system.

Referenced by solve_helper().

◆ C_pt

DoubleVector* oomph::AugmentedProblemGMRES::C_pt

private

Pointer to the row vector in the bordered system.

Referenced by augmented_matrix_multiply(), and solve_helper().

◆ Iteration_restart

bool oomph::AugmentedProblemGMRES::Iteration_restart

private

boolean indicating if iteration restarting is used

Referenced by AugmentedProblemGMRES(), disable_iteration_restart(), enable_iteration_restart(), iteration_restart(), and solve_helper().

◆ Iterations

unsigned oomph::AugmentedProblemGMRES::Iterations

private

Number of iterations taken.

Referenced by iterations(), and solve_helper().

◆ Matrix_can_be_deleted

bool oomph::AugmentedProblemGMRES::Matrix_can_be_deleted

private

Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.

Referenced by clean_up_memory(), and solve().

◆ Matrix_pt

CRDoubleMatrix* oomph::AugmentedProblemGMRES::Matrix_pt

private

Pointer to matrix.

Referenced by clean_up_memory(), resolve(), and solve().

◆ Preconditioner_LHS

bool oomph::AugmentedProblemGMRES::Preconditioner_LHS

private

boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false)

Referenced by AugmentedProblemGMRES(), set_preconditioner_LHS(), set_preconditioner_RHS(), solve_helper(), and update().

◆ Resolving

bool oomph::AugmentedProblemGMRES::Resolving

private

Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and preconditioner

Referenced by resolve(), solve(), and solve_helper().

◆ Restart

unsigned oomph::AugmentedProblemGMRES::Restart

private

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

◆ Rhs_pt

double* oomph::AugmentedProblemGMRES::Rhs_pt

private

Pointer to the last entry of the RHS vector in the bordered system.

◆ Schur_complement_scalar

double oomph::AugmentedProblemGMRES::Schur_complement_scalar

private

The scalar component of the Schur complement preconditioner.

Referenced by apply_schur_complement_preconditioner(), solve(), and solve_helper().

◆ X_pt

double* oomph::AugmentedProblemGMRES::X_pt

private

Pointer to the last entry of the LHS vector in the bordered system.

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

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ AugmentedProblemGMRES() [1/2]

◆ ~AugmentedProblemGMRES()

◆ AugmentedProblemGMRES() [2/2]

Member Function Documentation

◆ apply_plane_rotation()

◆ apply_schur_complement_preconditioner()

◆ augmented_matrix_multiply()

◆ clean_up_memory()

◆ disable_iteration_restart()

◆ disable_resolve()

◆ 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

◆ B_pt

◆ C_pt

◆ Iteration_restart

◆ Iterations

◆ Matrix_can_be_deleted

◆ Matrix_pt

◆ Preconditioner_LHS

◆ Resolving

◆ Restart

◆ Rhs_pt

◆ Schur_complement_scalar

◆ X_pt