#include <general_purpose_space_time_subsidiary_block_preconditioner.h>

Inheritance diagram for oomph::GMRESBlockPreconditioner:

Public Member Functions
	GMRESBlockPreconditioner ()
	Constructor (empty) More...

virtual	~GMRESBlockPreconditioner ()
	Destructor. More...

virtual void	clean_up_memory ()
	Clean up the memory (empty for now...) More...

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

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

void	setup ()
	Setup the preconditioner. More...

void	preconditioner_solve (const DoubleVector &r, DoubleVector &z)
	Apply preconditioner to r. More...

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

unsigned	iterations () const
	Handle to the number of iterations taken. More...

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

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

void	enable_doc_memory_usage ()
	Document the memory usage. More...

void	disable_doc_memory_usage ()
	Don't document the memory usage! More...

double	get_memory_usage_in_bytes ()
	Get the memory statistics. More...

SpaceTimeNavierStokesSubsidiaryPreconditioner *	navier_stokes_subsidiary_preconditioner_pt () const

	GMRESBlockPreconditioner ()
	Constructor (empty) More...

virtual	~GMRESBlockPreconditioner ()
	Destructor. More...

virtual void	clean_up_memory ()
	Clean up the memory (empty for now...) More...

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

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

void	setup ()
	Setup the preconditioner. More...

void	preconditioner_solve (const DoubleVector &r, DoubleVector &z)
	Apply preconditioner to r. More...

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

unsigned	iterations () const
	Handle to the number of iterations taken. More...

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

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

void	enable_doc_memory_statistics ()
	Document the memory usage. More...

void	disable_doc_memory_statistics ()
	Don't document the memory usage! More...

double	get_memory_usage_in_bytes ()
	Get the memory statistics. More...

SpaceTimeNavierStokesSubsidiaryPreconditioner *	navier_stokes_subsidiary_preconditioner_pt () const

void	setup (DoubleMatrixBase *matrix_pt)

void	setup (const Problem problem_pt, DoubleMatrixBase matrix_pt)

virtual void	setup ()=0

void	setup (DoubleMatrixBase *matrix_pt)

void	setup (const Problem problem_pt, DoubleMatrixBase matrix_pt)

virtual void	setup ()=0

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	disable_resolve ()

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

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

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 (const DoubleVector &rhs, DoubleVector &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)

Public Member Functions inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
	BlockPreconditioner ()
	Constructor. More...

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

virtual	~BlockPreconditioner ()
	Destructor. More...

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

CRDoubleMatrix *	matrix_pt () const

void	turn_on_recursive_debug_flag ()

void	turn_off_recursive_debug_flag ()

void	turn_on_debug_flag ()
	Toggles on the debug flag. More...

void	turn_off_debug_flag ()
	Toggles off the debug flag. More...

void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)

void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse, const Vector< Vector< unsigned >> &doftype_coarsen_map_coarse)

virtual void	block_setup ()

void	block_setup (const Vector< unsigned > &dof_to_block_map)

void	get_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const

CRDoubleMatrix	get_block (const unsigned &i, const unsigned &j, const bool &ignore_replacement_block=false) const

void	set_master_matrix_pt (CRDoubleMatrix *in_matrix_pt)
	Set the matrix_pt in the upper-most master preconditioner. More...

void	get_block_other_matrix (const unsigned &i, const unsigned &j, CRDoubleMatrix *in_matrix_pt, CRDoubleMatrix &output_matrix)

void	get_blocks (DenseMatrix< bool > &required_blocks, DenseMatrix< CRDoubleMatrix * > &block_matrix_pt) const

void	get_dof_level_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_block, const bool &ignore_replacement_block=false) const

void	get_dof_level_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix &output_block, const bool &ignore_replacement_block) const

CRDoubleMatrix	get_concatenated_block (const VectorMatrix< BlockSelector > &selected_block)

void	get_concatenated_block_vector (const Vector< unsigned > &block_vec_number, const DoubleVector &v, DoubleVector &b)

void	return_concatenated_block_vector (const Vector< unsigned > &block_vec_number, const DoubleVector &b, DoubleVector &v) const
	Takes concatenated block ordered vector, b, and copies its. More...

void	get_block_vectors (const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const

void	get_block_vectors (const DoubleVector &v, Vector< DoubleVector > &s) const

void	return_block_vectors (const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const

void	return_block_vectors (const Vector< DoubleVector > &s, DoubleVector &v) const

void	get_block_vector (const unsigned &n, const DoubleVector &v, DoubleVector &b) const

void	return_block_vector (const unsigned &n, const DoubleVector &b, DoubleVector &v) const

void	get_block_ordered_preconditioner_vector (const DoubleVector &v, DoubleVector &w)

void	return_block_ordered_preconditioner_vector (const DoubleVector &w, DoubleVector &v) const

unsigned	nblock_types () const
	Return the number of block types. More...

unsigned	ndof_types () const
	Return the total number of DOF types. More...

const Mesh *	mesh_pt (const unsigned &i) const

unsigned	nmesh () const

int	block_number (const unsigned &i_dof) const
	Return the block number corresponding to a global index i_dof. More...

int	index_in_block (const unsigned &i_dof) const

const LinearAlgebraDistribution *	block_distribution_pt (const unsigned &b) const
	Access function to the block distributions (const version). More...

LinearAlgebraDistribution *	block_distribution_pt (const unsigned b)
	Access function to the block distributions (non-const version). More...

LinearAlgebraDistribution *	dof_block_distribution_pt (const unsigned &b)
	Access function to the dof-level block distributions. More...

const LinearAlgebraDistribution *	master_distribution_pt () const

unsigned	ndof_types_in_mesh (const unsigned &i) const

bool	is_subsidiary_block_preconditioner () const

bool	is_master_block_preconditioner () const

void	set_block_output_to_files (const std::string &basefilename)

void	disable_block_output_to_files ()
	Turn off output of blocks (by clearing the basefilename string). More...

bool	block_output_on () const
	Test if output of blocks is on or not. More...

void	output_blocks_to_files (const std::string &basefilename, const unsigned &precision=8) const

void	post_block_matrix_assembly_partial_clear ()

BlockPreconditioner< CRDoubleMatrix > *	master_block_preconditioner_pt () const
	Access function to the master block preconditioner pt. More...

void	clear_block_preconditioner_base ()

void	document ()

Vector< Vector< unsigned > >	doftype_coarsen_map_fine () const

Vector< unsigned >	get_fine_grain_dof_types_in (const unsigned &i) const

unsigned	nfine_grain_dof_types_in (const unsigned &i) const

MapMatrix< unsigned, CRDoubleMatrix * >	replacement_dof_block_pt () const
	Access function to the replaced dof-level blocks. More...

void	setup_matrix_vector_product (MatrixVectorProduct matvec_prod_pt, CRDoubleMatrix block_pt, const Vector< unsigned > &block_col_indices)

void	setup_matrix_vector_product (MatrixVectorProduct matvec_prod_pt, CRDoubleMatrix block_pt, const unsigned &block_col_index)

void	internal_get_block_ordered_preconditioner_vector (const DoubleVector &v, DoubleVector &w) const

void	internal_return_block_ordered_preconditioner_vector (const DoubleVector &w, DoubleVector &v) const

unsigned	internal_nblock_types () const

unsigned	internal_ndof_types () const

void	internal_return_block_vector (const unsigned &n, const DoubleVector &b, DoubleVector &v) const

void	internal_get_block_vector (const unsigned &n, const DoubleVector &v, DoubleVector &b) const

void	internal_get_block_vectors (const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const

void	internal_get_block_vectors (const DoubleVector &v, Vector< DoubleVector > &s) const

void	internal_return_block_vectors (const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const

void	internal_return_block_vectors (const Vector< DoubleVector > &s, DoubleVector &v) const

void	internal_get_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_block) const

void	internal_get_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix &output_block) const

int	internal_block_number (const unsigned &i_dof) const

int	internal_index_in_block (const unsigned &i_dof) const

const LinearAlgebraDistribution *	internal_block_distribution_pt (const unsigned &b) const
	Access function to the internal block distributions. More...

void	insert_auxiliary_block_distribution (const Vector< unsigned > &block_vec_number, LinearAlgebraDistribution *dist_pt)

void	block_matrix_test (const unsigned &i, const unsigned &j, const CRDoubleMatrix *block_matrix_pt) const

int	get_index_of_value (const Vector< myType > &vec, const myType val, const bool sorted=false) const

Public Member Functions inherited from oomph::Preconditioner
	Preconditioner ()
	Constructor. More...

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

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

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

virtual void	preconditioner_solve_transpose (const DoubleVector &r, DoubleVector &z)

void	setup (DoubleMatrixBase *matrix_pt)

void	setup (const Problem problem_pt, DoubleMatrixBase matrix_pt)

void	enable_silent_preconditioner_setup ()
	Set up the block preconditioner quietly! More...

void	disable_silent_preconditioner_setup ()
	Be verbose in the block preconditioner setup. More...

virtual void	set_matrix_pt (DoubleMatrixBase *matrix_pt)
	Set the matrix pointer. More...

virtual const OomphCommunicator *	comm_pt () const
	Get function for comm pointer. More...

virtual void	set_comm_pt (const OomphCommunicator *const comm_pt)
	Set the communicator pointer. More...

double	setup_time () const
	Returns the time to setup the preconditioner. More...

virtual void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)

virtual void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse, const Vector< Vector< unsigned >> &doftype_coarsen_map_coarse)

Private Member Functions
void	update (const unsigned &k, const Vector< Vector< double >> &H, const Vector< double > &s, const Vector< DoubleVector > &v, const DoubleVector &block_x_with_size_of_full_x, 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)

void	update (const unsigned &k, const Vector< Vector< double >> &H, const Vector< double > &s, const Vector< DoubleVector > &v, const DoubleVector &block_x_with_size_of_full_x, 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
CRDoubleMatrix *	Matrix_pt
	Pointer to matrix. More...

SpaceTimeNavierStokesSubsidiaryPreconditioner *	Navier_stokes_subsidiary_preconditioner_pt
	Pointer to the preconditioner for the block matrix. More...

unsigned	Iterations
	Number of iterations taken. More...

bool	Compute_memory_statistics

double	Memory_usage_in_bytes

bool	Preconditioner_has_been_setup

bool	Preconditioner_LHS

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

Protected Member Functions inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
void	set_nmesh (const unsigned &n)

void	set_mesh (const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)

void	set_replacement_dof_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix *replacement_dof_block_pt)

bool	any_mesh_distributed () const

int	internal_dof_number (const unsigned &i_dof) const

unsigned	internal_index_in_dof (const unsigned &i_dof) const

unsigned	internal_block_dimension (const unsigned &b) const

unsigned	internal_dof_block_dimension (const unsigned &i) const

unsigned	master_nrow () const

unsigned	internal_master_dof_number (const unsigned &b) const

const LinearAlgebraDistribution *	internal_preconditioner_matrix_distribution_pt () const

const LinearAlgebraDistribution *	preconditioner_matrix_distribution_pt () const

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

Protected Attributes inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
MapMatrix< unsigned, CRDoubleMatrix * >	Replacement_dof_block_pt
	The replacement dof-level blocks. More...

Vector< LinearAlgebraDistribution * >	Block_distribution_pt
	The distribution for the blocks. More...

Vector< Vector< unsigned > >	Block_to_dof_map_coarse

Vector< Vector< unsigned > >	Block_to_dof_map_fine
	Mapping for the block types to the most fine grain dof types. More...

Vector< Vector< unsigned > >	Doftype_coarsen_map_coarse

Vector< Vector< unsigned > >	Doftype_coarsen_map_fine

Vector< LinearAlgebraDistribution * >	Internal_block_distribution_pt
	Storage for the default distribution for each internal block. More...

Vector< LinearAlgebraDistribution * >	Dof_block_distribution_pt

Vector< unsigned >	Allow_multiple_element_type_in_mesh

Vector< const Mesh * >	Mesh_pt

Vector< unsigned >	Ndof_types_in_mesh

unsigned	Internal_nblock_types

unsigned	Internal_ndof_types

Protected Attributes inherited from oomph::Preconditioner
bool	Silent_preconditioner_setup
	Boolean to indicate whether or not the build should be done silently. More...

std::ostream *	Stream_pt
	Pointer to the output stream – defaults to std::cout. More...

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

Detailed Description

The block preconditioner form of GMRES. This version extracts the blocks from the global systems and assembles the system by concatenating all the matrices together

Constructor & Destructor Documentation

◆ GMRESBlockPreconditioner() [1/4]

oomph::GMRESBlockPreconditioner::GMRESBlockPreconditioner ( )

inline

Constructor (empty)

       : BlockPreconditioner<CRDoubleMatrix>(),
         Matrix_pt(0),
         Navier_stokes_subsidiary_preconditioner_pt(0),
         Iterations(0),
         Preconditioner_has_been_setup(false),
         Preconditioner_LHS(false)
     {
       // By default, don't store the memory statistics of this preconditioner
       Compute_memory_statistics = false;
  
       // Initialise the value of Memory_usage_in_bytes
       Memory_usage_in_bytes = 0.0;
     } // End of GMRESBlockPreconditioner

References Compute_memory_statistics, and Memory_usage_in_bytes.

◆ ~GMRESBlockPreconditioner() [1/2]

virtual oomph::GMRESBlockPreconditioner::~GMRESBlockPreconditioner ( )

inlinevirtual

Destructor.

     {
       // Call the auxiliary clean up function
       this->clean_up_memory();
     } // End of ~GMRESBlockPreconditioner

References clean_up_memory().

◆ GMRESBlockPreconditioner() [2/4]

oomph::GMRESBlockPreconditioner::GMRESBlockPreconditioner ( const GMRESBlockPreconditioner & )

delete

Broken copy constructor.

◆ GMRESBlockPreconditioner() [3/4]

oomph::GMRESBlockPreconditioner::GMRESBlockPreconditioner ( )

inline

Constructor (empty)

       : BlockPreconditioner<CRDoubleMatrix>(),
         Matrix_pt(0),
         Navier_stokes_subsidiary_preconditioner_pt(0),
         Iterations(0),
         Preconditioner_has_been_setup(false),
         Preconditioner_LHS(false)
     {
       // By default, don't store the memory statistics of this preconditioner
       Compute_memory_statistics = false;
  
       // Initialise the value of Memory_usage_in_bytes
       Memory_usage_in_bytes = 0.0;
     } // End of GMRESBlockPreconditioner

References Compute_memory_statistics, and Memory_usage_in_bytes.

◆ ~GMRESBlockPreconditioner() [2/2]

virtual oomph::GMRESBlockPreconditioner::~GMRESBlockPreconditioner ( )

inlinevirtual

Destructor.

     {
       // Call the auxiliary clean up function
       this->clean_up_memory();
     } // End of ~GMRESBlockPreconditioner

References clean_up_memory().

◆ GMRESBlockPreconditioner() [4/4]

oomph::GMRESBlockPreconditioner::GMRESBlockPreconditioner ( const GMRESBlockPreconditioner & )

delete

Broken copy constructor.

Member Function Documentation

◆ apply_plane_rotation() [1/2]

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

Referenced by preconditioner_solve().

◆ apply_plane_rotation() [2/2]

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

◆ clean_up_memory() [1/2]

virtual void oomph::GMRESBlockPreconditioner::clean_up_memory ( )

inlinevirtual

Clean up the memory (empty for now...)

Reimplemented from oomph::LinearSolver.

     {
       // If the block preconditioner has been set up previously
       if (Preconditioner_has_been_setup)
       {
         // Delete the subsidiary preconditioner
         delete Navier_stokes_subsidiary_preconditioner_pt;
  
         // Make it a null pointer
         Navier_stokes_subsidiary_preconditioner_pt = 0;
  
         // Delete the matrix!
         delete Matrix_pt;
  
         // Make it a null pointer
         Matrix_pt = 0;
       }
     } // End of clean_up_memory

References Matrix_pt, Navier_stokes_subsidiary_preconditioner_pt, and Preconditioner_has_been_setup.

Referenced by setup(), and ~GMRESBlockPreconditioner().

◆ clean_up_memory() [2/2]

virtual void oomph::GMRESBlockPreconditioner::clean_up_memory ( )

inlinevirtual

Clean up the memory (empty for now...)

Reimplemented from oomph::LinearSolver.

     {
       // If the block preconditioner has been set up previously
       if (Preconditioner_has_been_setup)
       {
         // Delete the subsidiary preconditioner
         delete Navier_stokes_subsidiary_preconditioner_pt;
  
         // Make it a null pointer
         Navier_stokes_subsidiary_preconditioner_pt = 0;
  
         // Delete the matrix!
         delete Matrix_pt;
  
         // Make it a null pointer
         Matrix_pt = 0;
       }
     } // End of clean_up_memory

References Matrix_pt, Navier_stokes_subsidiary_preconditioner_pt, and Preconditioner_has_been_setup.

◆ disable_doc_memory_statistics()

void oomph::GMRESBlockPreconditioner::disable_doc_memory_statistics ( )

inline

Don't document the memory usage!

Set the appropriate flag to false

     {
       Compute_memory_statistics = false;
     } // End of disable_doc_memory_statistics

References Compute_memory_statistics.

◆ disable_doc_memory_usage()

void oomph::GMRESBlockPreconditioner::disable_doc_memory_usage ( )

inline

Don't document the memory usage!

Set the appropriate flag to false

     {
       Compute_memory_statistics = false;
     } // End of disable_doc_memory_usage

References Compute_memory_statistics.

◆ enable_doc_memory_statistics()

void oomph::GMRESBlockPreconditioner::enable_doc_memory_statistics ( )

inline

Document the memory usage.

Set the appropriate flag to true

     {
       Compute_memory_statistics = true;
     } // End of enable_doc_memory_statistics

References Compute_memory_statistics.

◆ enable_doc_memory_usage()

void oomph::GMRESBlockPreconditioner::enable_doc_memory_usage ( )

inline

Document the memory usage.

Set the appropriate flag to true

     {
       Compute_memory_statistics = true;
     } // End of enable_doc_memory_usage

References Compute_memory_statistics.

◆ generate_plane_rotation() [1/2]

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

◆ generate_plane_rotation() [2/2]

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

◆ get_memory_usage_in_bytes() [1/2]

double oomph::GMRESBlockPreconditioner::get_memory_usage_in_bytes ( )

inline

Get the memory statistics.

     {
       // Has the preconditioner even been set up yet?
       if (Preconditioner_has_been_setup)
       {
         // Were we meant to compute the statistics?
         if (Compute_memory_statistics)
         {
           // Return the appropriate variable value
           return Memory_usage_in_bytes;
         }
         else
         {
           // Allocate storage for an output stream
           std::ostringstream warning_message_stream;
  
           // Create a warning message
           warning_message_stream
             << "The memory statistics have not been calculated "
             << "so I'm returning\nthe value zero." << std::endl;
  
           // Give the user a warning
           OomphLibWarning(warning_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
  
           // Return the value zero
           return 0.0;
         }
       }
       // If the preconditioner hasn't been set up yet
       else
       {
         // Allocate storage for an output stream
         std::ostringstream warning_message_stream;
  
         // Create a warning message
         warning_message_stream
           << "The preconditioner hasn't even been set up yet "
           << "so I'm returning\nthe value zero." << std::endl;
  
         // Give the user a warning
         OomphLibWarning(warning_message_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
  
         // Return the value zero
         return 0.0;
       } // if (Preconditioner_has_been_setup)
     } // End of get_memory_usage_in_bytes

References Compute_memory_statistics, Memory_usage_in_bytes, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Preconditioner_has_been_setup.

Referenced by oomph::BandedBlockTriangularPreconditioner< MATRIX >::setup().

◆ get_memory_usage_in_bytes() [2/2]

double oomph::GMRESBlockPreconditioner::get_memory_usage_in_bytes ( )

inline

Get the memory statistics.

     {
       // Has the preconditioner even been set up yet?
       if (Preconditioner_has_been_setup)
       {
         // Were we meant to compute the statistics?
         if (Compute_memory_statistics)
         {
           // Return the appropriate variable value
           return Memory_usage_in_bytes;
         }
         else
         {
           // Allocate storage for an output stream
           std::ostringstream warning_message_stream;
  
           // Create a warning message
           warning_message_stream
             << "The memory statistics have not been calculated "
             << "so I'm returning\nthe value zero." << std::endl;
  
           // Give the user a warning
           OomphLibWarning(warning_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
  
           // Return the value zero
           return 0.0;
         }
       }
       // If the preconditioner hasn't been set up yet
       else
       {
         // Allocate storage for an output stream
         std::ostringstream warning_message_stream;
  
         // Create a warning message
         warning_message_stream
           << "The preconditioner hasn't even been set up yet "
           << "so I'm returning\nthe value zero." << std::endl;
  
         // Give the user a warning
         OomphLibWarning(warning_message_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
  
         // Return the value zero
         return 0.0;
       } // if (Preconditioner_has_been_setup)
     } // End of get_memory_usage_in_bytes

References Compute_memory_statistics, Memory_usage_in_bytes, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Preconditioner_has_been_setup.

◆ iterations() [1/2]

unsigned oomph::GMRESBlockPreconditioner::iterations ( ) const

inlinevirtual

Handle to the number of iterations taken.

Implements oomph::IterativeLinearSolver.

     {
       // Return the number of iterations
       return Iterations;
     } // End of iterations

References Iterations.

◆ iterations() [2/2]

unsigned oomph::GMRESBlockPreconditioner::iterations ( ) const

inlinevirtual

Handle to the number of iterations taken.

Implements oomph::IterativeLinearSolver.

     {
       // Return the number of iterations
       return Iterations;
     } // End of iterations

References Iterations.

◆ navier_stokes_subsidiary_preconditioner_pt() [1/2]

SpaceTimeNavierStokesSubsidiaryPreconditioner* oomph::GMRESBlockPreconditioner::navier_stokes_subsidiary_preconditioner_pt ( ) const

inline

Handle to the Navier-Stokes subsidiary block preconditioner DRAIG: Make sure the desired const-ness is correct later...

     {
       // Return a pointer to the appropriate member data
       return Navier_stokes_subsidiary_preconditioner_pt;
     } // End of navier_stokes_subsidiary_preconditioner_pt

References Navier_stokes_subsidiary_preconditioner_pt.

Referenced by oomph::BandedBlockTriangularPreconditioner< MATRIX >::setup().

◆ navier_stokes_subsidiary_preconditioner_pt() [2/2]

SpaceTimeNavierStokesSubsidiaryPreconditioner* oomph::GMRESBlockPreconditioner::navier_stokes_subsidiary_preconditioner_pt ( ) const

inline

Handle to the Navier-Stokes subsidiary block preconditioner DRAIG: Make sure the desired const-ness is correct later...

     {
       // Return a pointer to the appropriate member data
       return Navier_stokes_subsidiary_preconditioner_pt;
     } // End of navier_stokes_subsidiary_preconditioner_pt

References Navier_stokes_subsidiary_preconditioner_pt.

◆ operator=() [1/2]

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

delete

Broken assignment operator.

◆ operator=() [2/2]

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

delete

Broken assignment operator.

◆ preconditioner_solve() [1/2]

void oomph::GMRESBlockPreconditioner::preconditioner_solve	(	const DoubleVector &	r,
		DoubleVector &	z
	)

virtual

Apply preconditioner to r.

Preconditioner solve for the GMRES block preconditioner.

Implements oomph::Preconditioner.

   {
     // Get the number of block types
     unsigned n_block_types = this->nblock_types();
  
     // If there is more than one block type (there shouldn't be!)
     if (n_block_types != 1)
     {
       // Throw an error
       throw OomphLibError("Can only deal with one dof type!",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
     // Allocate space for the solution blocks
     Vector<DoubleVector> block_z(n_block_types);
  
     // Allocate space for the residual blocks
     Vector<DoubleVector> block_r(n_block_types);
  
     // Get the reordered RHS vector
     this->get_block_vectors(rhs, block_r);
  
     // Get the reordered solution vector
     this->get_block_vectors(solution, block_z);
  
     // Initialise the entries in the solution vector to zero
     block_z[0].initialise(0.0);
  
     // Get number of dofs
     unsigned n_dof = block_r[0].nrow();
  
     // Time solver
     double t_start = TimingHelpers::timer();
  
     // Relative residual
     double resid;
  
     // Iteration counter
     unsigned iter = 1;
  
     // Initialise vectors
     Vector<double> s(n_dof + 1, 0);
     Vector<double> cs(n_dof + 1);
     Vector<double> sn(n_dof + 1);
     DoubleVector w(block_r[0].distribution_pt(), 0.0);
  
     // Storage for the time taken for all preconditioner applications
     double t_prec_application_total = 0.0;
  
     // Storage for the RHS vector
     DoubleVector r(block_r[0].distribution_pt(), 0.0);
  
     // If we're using LHS preconditioning then solve Mr=b-Jx for r (assumes x=0)
     if (Preconditioner_LHS)
     {
       // Initialise timer
       double t_prec_application_start = TimingHelpers::timer();
  
       // This is pretty inefficient but there is no other way to do this with
       // block subsidiary preconditioners as far as I can tell because they
       // demand to have the full r and z vectors...
       DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(), 0.0);
       DoubleVector r_updated(rhs.distribution_pt(), 0.0);
  
       // Reorder the entries in block_r[0] and put them back into the
       // global-sized vector. The subsidiary preconditioner should only
       // operate on the dofs that this preconditioner operates on so
       // don't worry about all the other entries
       this->return_block_vector(0, block_r[0], r_updated);
  
       // Solve on the global-sized vectors
       Navier_stokes_subsidiary_preconditioner_pt->preconditioner_solve(
         r_updated, block_z_with_size_of_full_z);
  
       // Now grab the part of the solution associated with this preconditioner
       this->get_block_vector(0, block_z_with_size_of_full_z, block_z[0]);
  
       // 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, do nothing to the actual RHS vector
     else
     {
       // Copy the entries into r
       r = block_r[0];
     }
  
     // Calculate the initial residual norm (assumes x=0 initially)
     double normb = r.norm();
  
     // Set beta (the initial residual)
     double beta = normb;
  
     // Compute initial relative residual
     if (normb == 0.0)
     {
       // To avoid dividing by zero, set normb to 1
       normb = 1;
     }
  
     // Now calculate the relative residual norm
     resid = beta / normb;
  
     // If required, document convergence history to screen or file (if
     // stream is open)
     if (Doc_convergence_history)
     {
       // If a file has not been provided to output to
       if (!Output_file_stream.is_open())
       {
         // Output the initial relative residual norm to screen
         oomph_info << 0 << " " << resid << std::endl;
       }
       // If there is a file provided to output to
       else
       {
         // Output the initial relative residual norm to file
         Output_file_stream << 0 << " " << resid << std::endl;
       }
     } // if (Doc_convergence_history)
  
     // If GMRES converges immediately
     if (resid <= Tolerance)
     {
       // If the user wishes GMRES to output information about the solve to
       // screen
       if (Doc_time)
       {
         // Tell the user
         oomph_info << "GMRES block preconditioner converged immediately. "
                    << "Normalised residual norm: " << resid << std::endl;
       }
  
       // Now reorder the values in block_z[0] and put them back into the
       // global solution vector
       this->return_block_vector(0, block_z[0], solution);
  
       // Doc time for solver
       double t_end = TimingHelpers::timer();
       Solution_time = t_end - t_start;
  
       // If the user wishes GMRES to output information about the solve to
       // screen
       if (Doc_time)
       {
         // Tell the user
         oomph_info << "Time for all preconditioner applications [sec]: "
                    << t_prec_application_total
                    << "\nTotal time for solve with GMRES block preconditioner "
                    << "[sec]: " << Solution_time << std::endl;
       }
  
       // We're finished; end here!
       return;
     } // if (resid<=Tolerance)
  
     // Initialise vector of orthogonal basis vectors, v
     Vector<DoubleVector> v(n_dof + 1);
  
     // Initialise the upper Hessenberg matrix H.
     // NOTE: For implementation purposes the upper hessenberg matrix indices
     // are swapped so the matrix is effectively transposed
     Vector<Vector<double>> H(n_dof + 1);
  
     // Loop as many times as we're allowed
     while (iter <= Max_iter)
     {
       // Copy the entries of the residual vector, r, into v[0]
       v[0] = r;
  
       // Now update the zeroth basis vector, v[0], to have values r/beta
       v[0] /= beta;
  
       // Storage for ...?
       s[0] = beta;
  
       // Inner iteration counter for restarted version (we don't actually use
       // a restart in this implementation of GMRES...)
       unsigned iter_restart;
  
       // Perform iterations
       for (iter_restart = 0; iter_restart < n_dof && iter <= Max_iter;
            iter_restart++, iter++)
       {
         // Resize next column of the upper Hessenberg matrix
         H[iter_restart].resize(iter_restart + 2);
  
         // Solve for w inside the scope of these braces
         {
           // Temporary vector
           DoubleVector temp(block_r[0].distribution_pt(), 0.0);
  
           // If we're using LHS preconditioning solve Mw=Jv[i] for w
           if (Preconditioner_LHS)
           {
             // Calculate temp=Jv[i]
             Matrix_pt->multiply(v[iter_restart], temp);
  
             // Get the current time
             double t_prec_application_start = TimingHelpers::timer();
  
             // This is pretty inefficient but there is no other way to do this
             // with block sub preconditioners as far as I can tell because they
             // demand to have the full r and z vectors...
             DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(),
                                                      0.0);
             DoubleVector r_updated(rhs.distribution_pt(), 0.0);
  
             // Put the values in temp into the global-sized vector
             this->return_block_vector(0, temp, r_updated);
  
             // Apply the NSSP to the global-sized vectors (this will extract
             // the appropriate sub-vectors and solve on the reduced system)
             Navier_stokes_subsidiary_preconditioner_pt->preconditioner_solve(
               r_updated, block_z_with_size_of_full_z);
  
             // Now grab the part of the solution associated with this
             // preconditioner
             this->get_block_vector(0, block_z_with_size_of_full_z, 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);
           }
           // If we're using RHS preconditioning solve
           else
           {
             // Initialise timer
             double t_prec_application_start = TimingHelpers::timer();
  
             // This is pretty inefficient but there is no other way to do this
             // with block sub preconditioners as far as I can tell because they
             // demand to have the full r and z vectors...
             DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt());
             DoubleVector r_updated(rhs.distribution_pt());
  
             // Put the values in v[iter_restart] back into the global-sized
             // vector
             this->return_block_vector(0, v[iter_restart], r_updated);
  
             // Solve on the global-sized vectors
             Navier_stokes_subsidiary_preconditioner_pt->preconditioner_solve(
               r_updated, block_z_with_size_of_full_z);
  
             // Now grab the part of the solution associated with this
             // preconditioner
             this->get_block_vector(0, block_z_with_size_of_full_z, 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
             Matrix_pt->multiply(temp, w);
           }
         } // End of solve for w
  
         // Get a pointer to the values in w
         double* w_pt = w.values_pt();
  
         // Loop over the columns of the (transposed) Hessenberg matrix
         for (unsigned k = 0; k <= iter_restart; k++)
         {
           // Initialise the entry in the upper Hessenberg matrix
           H[iter_restart][k] = 0.0;
  
           // Get the pointer to the values of the k-th basis vector
           double* vk_pt = v[k].values_pt();
  
           // Update H with the dot product of w and v[k]
           H[iter_restart][k] += w.dot(v[k]);
  
           // Loop over the entries in w and v[k] again
           for (unsigned i = 0; i < n_dof; i++)
           {
             // Now update the values in w
             w_pt[i] -= H[iter_restart][k] * vk_pt[i];
           }
         } // for (unsigned k=0;k<=iter_restart;k++)
  
         // Calculate the (iter_restart+1,iter_restart)-th entry in the upper
         // Hessenberg matrix (remember, H is actually the transpose of the
         // proper upper Hessenberg matrix)
         H[iter_restart][iter_restart + 1] = w.norm();
  
         // Build the (iter_restart+1)-th basis vector by copying w
         v[iter_restart + 1] = w;
  
         // Now rescale the basis vector
         v[iter_restart + 1] /= H[iter_restart][iter_restart + 1];
  
         // Loop over the columns of the transposed upper Hessenberg matrix
         for (unsigned k = 0; k < iter_restart; k++)
         {
           // Apply a Givens rotation to entries of H
           apply_plane_rotation(
             H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
         }
  
         // Now generate the cs and sn entries for a plane rotation
         generate_plane_rotation(H[iter_restart][iter_restart],
                                 H[iter_restart][iter_restart + 1],
                                 cs[iter_restart],
                                 sn[iter_restart]);
  
         // Use the newly generated entries in cs and sn to apply a Givens
         // rotation to those same entries
         apply_plane_rotation(H[iter_restart][iter_restart],
                              H[iter_restart][iter_restart + 1],
                              cs[iter_restart],
                              sn[iter_restart]);
  
         // Now apply the Givens rotation to the entries in s
         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 the relative residual norm
         resid = beta / normb;
  
         // If required, document the convergence history to screen or file
         if (Doc_convergence_history)
         {
           // If an output file has not been provided
           if (!Output_file_stream.is_open())
           {
             // Output the convergence history to screen
             oomph_info << iter << " " << resid << std::endl;
           }
           // An output file has been provided so output to that
           else
           {
             // Output the convergence history to file
             Output_file_stream << iter << " " << resid << std::endl;
           }
         } // if (Doc_convergence_history)
  
         // If the required tolerance was met
         if (resid < Tolerance)
         {
           // Storage for the global-sized solution vector. Strictly speaking
           // we could actually just use the vector, solution but this can be
           // done after we know we've implemented everything correctly...
           DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(), 0.0);
  
           // Take the current solution, reorder the entries appropriately
           // and stick them into the global sized vector
           this->return_block_vector(0, block_z[0], block_z_with_size_of_full_z);
  
           // This will update block_z[0] and won't touch the global-size vector
           // block_x_with_size_of_full_x. We're in charge of reordering the
           // entries and putting it in solution
           update(
             iter_restart, H, s, v, block_z_with_size_of_full_z, block_z[0]);
  
           // Now go into block_z[0], reorder the entries in the appropriate
           // manner and put them into the vector, solution
           this->return_block_vector(0, block_z[0], solution);
  
           // If the user wishes GMRES to document the convergence
           if (Doc_time)
           {
             // Document the convergence to screen
             oomph_info
               << "\nGMRES block preconditioner converged (1). Normalised "
               << "residual norm: " << resid
               << "\nNumber of iterations to convergence: " << iter << "\n"
               << std::endl;
           }
  
           // Get the current time
           double t_end = TimingHelpers::timer();
  
           // Calculate the time difference, i.e. the time taken for the whole
           // solve
           Solution_time = t_end - t_start;
  
           // Storage the iteration count
           Iterations = iter;
  
           // If the user wishes GMRES to document the convergence
           if (Doc_time)
           {
             // Document the convergence to screen
             oomph_info
               << "Time for all preconditioner applications [sec]: "
               << t_prec_application_total
               << "\nTotal time for solve with GMRES block preconditioner "
               << "[sec]: " << Solution_time << std::endl;
           }
  
           // We're done now; finish here
           return;
         } // if (resid<Tolerance)
       } // for (iter_restart=0;iter_restart<n_dof&&iter<=Max_iter;...
  
       // Update
       if (iter_restart > 0)
       {
         // Storage for the global-sized solution vector
         DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(), 0.0);
  
         // Take the current solution, reorder the entries appropriately
         // and stick them into the global sized vector
         this->return_block_vector(0, block_z[0], block_z_with_size_of_full_z);
  
         // Update the (reordered block) solution vector
         update(
           (iter_restart - 1), H, s, v, block_z_with_size_of_full_z, block_z[0]);
  
         // Now go into block_z[0], reorder the entries in the appropriate manner
         // and put them into the vector, solution
         this->return_block_vector(0, block_z[0], solution);
       }
  
       // Solve Mr=(b-Jx) for r
       {
         // Temporary vector to store b-Jx
         DoubleVector temp(block_r[0].distribution_pt(), 0.0);
  
         // Calculate temp=b-Jx
         Matrix_pt->residual(block_z[0], block_r[0], temp);
  
         // If we're using LHS preconditioning
         if (Preconditioner_LHS)
         {
           // Initialise timer
           double t_prec_application_start = TimingHelpers::timer();
  
           // This is pretty inefficient but there is no other way to do this
           // with block sub preconditioners as far as I can tell because they
           // demand to have the full r and z vectors...
           DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt());
           DoubleVector r_updated(rhs.distribution_pt());
  
           // Reorder the entries of temp and put them into the global-sized
           // vector
           this->return_block_vector(0, temp, r_updated);
  
           // Solve on the global-sized vectors
           Navier_stokes_subsidiary_preconditioner_pt->preconditioner_solve(
             r_updated, block_z_with_size_of_full_z);
  
           // Now grab the part of the solution associated with this
           // preconditioner
           this->get_block_vector(0, block_z_with_size_of_full_z, 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);
         }
       } // End of solve Mr=(b-Jx) for r
  
       // Compute the current residual norm
       beta = r.norm();
  
       // if relative residual within tolerance
       resid = beta / normb;
       if (resid < Tolerance)
       {
         // Now store the result in the global solution vector
         this->return_block_vector(0, block_z[0], solution);
  
         if (Doc_time)
         {
           oomph_info
             << "\nGMRES block preconditioner converged (2). 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 (Doc_time)
         {
           oomph_info
             << "Time for all preconditioner applications [sec]: "
             << t_prec_application_total
             << "\nTotal time for solve with GMRES block preconditioner "
             << "[sec]: " << Solution_time << std::endl;
         }
         return;
       }
     }
  
  
     // otherwise GMRES failed convergence
     oomph_info << "\nGMRES block preconditioner 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::ostringstream error_message_stream;
       error_message_stream << "Solver failed to converge and you requested an "
                            << "error on convergence failures.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_EXCEPTION_LOCATION,
                           OOMPH_CURRENT_FUNCTION);
     }
   } // End of solve_helper

References apply_plane_rotation(), beta, oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::IterativeLinearSolver::Doc_convergence_history, oomph::LinearSolver::Doc_time, boost::multiprecision::fabs(), generate_plane_rotation(), oomph::BlockPreconditioner< CRDoubleMatrix >::get_block_vector(), oomph::BlockPreconditioner< CRDoubleMatrix >::get_block_vectors(), H, i, oomph::Vector< _Tp >::initialise(), Iterations, k, Matrix_pt, oomph::IterativeLinearSolver::Max_iter, oomph::CRDoubleMatrix::multiply(), Navier_stokes_subsidiary_preconditioner_pt, oomph::BlockPreconditioner< CRDoubleMatrix >::nblock_types(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::IterativeLinearSolver::Output_file_stream, Preconditioner_LHS, oomph::SpaceTimeNavierStokesSubsidiaryPreconditioner::preconditioner_solve(), UniformPSDSelfTest::r, oomph::DoubleMatrixBase::residual(), oomph::BlockPreconditioner< CRDoubleMatrix >::return_block_vector(), s, BiharmonicTestFunctions1::solution(), oomph::IterativeLinearSolver::Solution_time, oomph::IterativeLinearSolver::Throw_error_after_max_iter, oomph::TimingHelpers::timer(), oomph::IterativeLinearSolver::Tolerance, update(), v, and w.

◆ preconditioner_solve() [2/2]

void oomph::GMRESBlockPreconditioner::preconditioner_solve	(	const DoubleVector &	r,
		DoubleVector &	z
	)

virtual

Apply preconditioner to r.

Implements oomph::Preconditioner.

◆ set_preconditioner_LHS() [1/2]

void oomph::GMRESBlockPreconditioner::set_preconditioner_LHS ( )

inline

Set left preconditioning (the default)

     {
       Preconditioner_LHS = true;
     }

References Preconditioner_LHS.

◆ set_preconditioner_LHS() [2/2]

void oomph::GMRESBlockPreconditioner::set_preconditioner_LHS ( )

inline

Set left preconditioning (the default)

     {
       Preconditioner_LHS = true;
     }

References Preconditioner_LHS.

◆ set_preconditioner_RHS() [1/2]

void oomph::GMRESBlockPreconditioner::set_preconditioner_RHS ( )

inline

Enable right preconditioning.

     {
       Preconditioner_LHS = false;
     }

References Preconditioner_LHS.

◆ set_preconditioner_RHS() [2/2]

void oomph::GMRESBlockPreconditioner::set_preconditioner_RHS ( )

inline

Enable right preconditioning.

     {
       Preconditioner_LHS = false;
     }

References Preconditioner_LHS.

◆ setup() [1/8]

void oomph::GMRESBlockPreconditioner::setup ( )

virtual

Setup the preconditioner.

Setup for the GMRES block preconditioner.

Implements oomph::Preconditioner.

   {
     // Clean up any previously created data
     clean_up_memory();
  
     // If we're meant to build silently
     if (this->Silent_preconditioner_setup == true)
     {
       // Store the output stream pointer
       this->Stream_pt = oomph_info.stream_pt();
  
       // Now set the oomph_info stream pointer to the null stream to
       // disable all possible output
       oomph_info.stream_pt() = &oomph_nullstream;
     }
  
     // Start the timer for the block setup
     double t_block_setup_start = TimingHelpers::timer();
  
     // Set up block look up schemes (done automatically in the
     // BlockPreconditioner base class, based on the information
     // provided in the block-preconditionable elements in the problem)
  
     // This preconditioner has two types of block:
     //        Type 0: Velocity - Corresponding to DOF type 0 & 1
     //        Type 1: Pressure - Corresponding to DOF type 2
     unsigned n_dof_types = 0;
  
     // If this has been set up as a subsidiary preconditioner
     if (this->is_subsidiary_block_preconditioner())
     {
       // Get the number of dof types (will be told by the master preconditioner)
       n_dof_types = this->ndof_types();
  
       // Output the number of dof types
       oomph_info << "Number of dof types: " << n_dof_types << std::endl;
  
       // Make sure there's only 3 dof types for now!
       if (n_dof_types != 3)
       {
         // Allocate space for an error message
         std::ostringstream error_message;
  
         // Create the error message
         error_message << "Should only be 3 dof types! You have " << n_dof_types
                       << " dof types!" << std::endl;
  
         // Throw an error
         throw OomphLibError(error_message.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
     // If it's being used as a master preconditioner
     else
     {
       // Throw an error
       throw OomphLibError("Currently only used as a subsidiary preconditioner!",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     } // if (this->is_subsidiary_block_preconditioner())
  
     // Aggregrate everything together since GMRES itself isn't going to take
     // advantage of the block structure...
     Vector<unsigned> dof_to_block_map(n_dof_types, 0);
  
     // Call the generic block setup function
     this->block_setup(dof_to_block_map);
  
     // Stop the timer
     double t_block_setup_finish = TimingHelpers::timer();
  
     // Calculate the time difference
     double block_setup_time = t_block_setup_finish - t_block_setup_start;
  
     // Document the time taken
     oomph_info << "Time for block_setup(...) [sec]: " << block_setup_time
                << std::endl;
  
     // Get the number of block types
     unsigned n_block_type = this->nblock_types();
  
     // How many block types are there? Should only be 2...
     oomph_info << "Number of block types: " << n_block_type << std::endl;
  
     // Space for the concatenated block matrix
     Matrix_pt = new CRDoubleMatrix;
  
     // Get the (one and only) block to be used by this block preconditioner
     this->get_block(0, 0, *Matrix_pt);
  
     // Create a new instance of the NS subsidiary preconditioner
     Navier_stokes_subsidiary_preconditioner_pt =
       new SpaceTimeNavierStokesSubsidiaryPreconditioner;
  
     // Do we need to doc. the memory statistics?
     if (Compute_memory_statistics)
     {
       // How many rows are there in this block?
       unsigned n_row = Matrix_pt->nrow();
  
       // How many nonzeros are there in this block?
       unsigned n_nnz = Matrix_pt->nnz();
  
       // Compute the memory usage. The only memory overhead here (for the
       // setup phase) is storing the system matrix
       Memory_usage_in_bytes = ((2 * ((n_row + 1) * sizeof(int))) +
                                (n_nnz * (sizeof(int) + sizeof(double))));
  
       // We might as well document the memory statistics of the Navier-Stokes
       // subsidiary block preconditioner while we're at it...
       Navier_stokes_subsidiary_preconditioner_pt->enable_doc_memory_usage();
     }
  
     // Create a map to declare which of the 3 dof types in the GMRES
     // subsidiary preconditioner correspond to the dof types in the NS
     // subsidiary block preconditioner
     Vector<unsigned> dof_map(n_dof_types);
  
     // Loop over the dof types associated with the GMRES subsidiary block
     // preconditioner
     for (unsigned i = 0; i < n_dof_types; i++)
     {
       // There is, essentially, an identity mapping between the GMRES subsidiary
       // block preconditioner and the NSSP w.r.t. the dof types
       dof_map[i] = i;
     }
  
     // Now turn it into an "proper" subsidiary preconditioner
     Navier_stokes_subsidiary_preconditioner_pt
       ->turn_into_subsidiary_block_preconditioner(this, dof_map);
  
     // Setup: The NS subsidiary preconditioner is itself a block preconditioner
     // so we need to pass it a pointer to full-size matrix!
     Navier_stokes_subsidiary_preconditioner_pt->setup(this->matrix_pt());
  
     // Remember that the preconditioner has been setup so the stored
     // information can be wiped when we come here next...
     Preconditioner_has_been_setup = true;
  
     // If we're meant to build silently, reassign the oomph stream pointer
     if (this->Silent_preconditioner_setup == true)
     {
       // Store the output stream pointer
       oomph_info.stream_pt() = this->Stream_pt;
  
       // Reset our own stream pointer
       this->Stream_pt = 0;
     }
   } // End of setup

References oomph::BlockPreconditioner< CRDoubleMatrix >::block_setup(), clean_up_memory(), Compute_memory_statistics, oomph::SpaceTimeNavierStokesSubsidiaryPreconditioner::enable_doc_memory_usage(), oomph::BlockPreconditioner< CRDoubleMatrix >::get_block(), i, int(), oomph::BlockPreconditioner< CRDoubleMatrix >::is_subsidiary_block_preconditioner(), oomph::BlockPreconditioner< CRDoubleMatrix >::matrix_pt(), Matrix_pt, Memory_usage_in_bytes, Navier_stokes_subsidiary_preconditioner_pt, oomph::BlockPreconditioner< CRDoubleMatrix >::nblock_types(), oomph::BlockPreconditioner< CRDoubleMatrix >::ndof_types(), oomph::CRDoubleMatrix::nnz(), oomph::CRDoubleMatrix::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::oomph_nullstream, Preconditioner_has_been_setup, oomph::SpaceTimeNavierStokesSubsidiaryPreconditioner::setup(), oomph::Preconditioner::Silent_preconditioner_setup, oomph::OomphInfo::stream_pt(), oomph::Preconditioner::Stream_pt, oomph::TimingHelpers::timer(), and oomph::BlockPreconditioner< MATRIX >::turn_into_subsidiary_block_preconditioner().

◆ setup() [2/8]

void oomph::GMRESBlockPreconditioner::setup ( )

virtual

Setup the preconditioner.

Implements oomph::Preconditioner.