The GMRES method for the Helmholtz solver. More...

#include <complex_smoother.h>

Inheritance diagram for oomph::HelmholtzGMRESMG< MATRIX >:

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

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

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

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

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

void	preconditioner_solve (const DoubleVector &r, DoubleVector &z)

void	setup ()

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

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)

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

virtual	~BlockPreconditioner ()
	Destructor. More...

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

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

MATRIX *	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< MATRIX > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)

void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< MATRIX > *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, MATRIX &output_matrix, const bool &ignore_replacement_block=false) const

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

void	set_master_matrix_pt (MATRIX *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, MATRIX *in_matrix_pt, MATRIX &output_matrix)

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

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

MATRIX	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< MATRIX > *	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, MATRIX &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 MATRIX *block_matrix_pt) const

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

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

void	get_dof_level_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix &output_block, const bool &ignore_replacement_block) 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)

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

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

void	complex_matrix_multiplication (Vector< CRDoubleMatrix * > const matrices_pt, const Vector< DoubleVector > &x, Vector< DoubleVector > &soln)

void	update (const unsigned &k, const Vector< Vector< std::complex< double >>> &hessenberg, const Vector< std::complex< double >> &s, const Vector< Vector< DoubleVector >> &v, Vector< DoubleVector > &x)
	Helper function to update the result vector. More...

void	generate_plane_rotation (std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)

void	apply_plane_rotation (std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)

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

Protected Member Functions inherited from oomph::BlockPreconditioner< MATRIX >
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
unsigned	Iterations
	Number of iterations taken. More...

Vector< CRDoubleMatrix * >	Matrices_storage_pt
	Vector of pointers to the real and imaginary part of the system matrix. More...

bool	Resolving

bool	Matrix_can_be_deleted

bool	Preconditioner_LHS

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< MATRIX >
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...

Additional Inherited Members
Static Protected Attributes inherited from oomph::IterativeLinearSolver
static IdentityPreconditioner	Default_preconditioner

Detailed Description

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

The GMRES method for the Helmholtz solver.

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

Constructor & Destructor Documentation

◆ HelmholtzGMRESMG() [1/2]

template<typename MATRIX >

oomph::HelmholtzGMRESMG< MATRIX >::HelmholtzGMRESMG ( )

inline

Constructor.

       : BlockPreconditioner<CRDoubleMatrix>(),
         Iterations(0),
         Resolving(false),
         Matrix_can_be_deleted(true)
     {
       Preconditioner_LHS = true;
     }

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

◆ ~HelmholtzGMRESMG()

template<typename MATRIX >

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

inlinevirtual

Destructor (cleanup storage)

     {
       clean_up_memory();
     }

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

◆ HelmholtzGMRESMG() [2/2]

template<typename MATRIX >

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

delete

Broken copy constructor.

Member Function Documentation

◆ apply_plane_rotation()

template<typename MATRIX >

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

inlineprotected

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

\[ \begin{bmatrix} dx \\ dy \end{bmatrix} \leftarrow \begin{bmatrix} \overline{\cos\theta} & \overline{\sin\theta} \\ -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \\ dy \end{bmatrix}. \]

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

     {
       // Calculate the value of dx but don't update it yet
       std::complex<double> temp = std::conj(cs) * dx + std::conj(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

References conj().

◆ clean_up_memory()

template<typename MATRIX >

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

inlineprotectedvirtual

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

Reimplemented from oomph::LinearSolver.

     {
       // If the matrix storage has been resized
       if (Matrices_storage_pt.size() > 0)
       {
         // If the real matrix pointer isn't null AND we're allowed to delete
         // the matrix which is only when we create the matrix ourselves
         if ((Matrices_storage_pt[0] != 0) && (Matrix_can_be_deleted))
         {
           // Delete the matrix
           delete Matrices_storage_pt[0];
  
           // Assign the associated pointer the value NULL
           Matrices_storage_pt[0] = 0;
         }
  
         // If the real matrix pointer isn't null AND we're allowed to delete
         // the matrix which is only when we create the matrix ourselves
         if ((Matrices_storage_pt[1] != 0) && (Matrix_can_be_deleted))
         {
           // Delete the matrix
           delete Matrices_storage_pt[1];
  
           // Assign the associated pointer the value NULL
           Matrices_storage_pt[1] = 0;
         }
       }
     } // End of clean_up_memory

References oomph::HelmholtzGMRESMG< MATRIX >::Matrices_storage_pt, and oomph::HelmholtzGMRESMG< MATRIX >::Matrix_can_be_deleted.

Referenced by oomph::HelmholtzGMRESMG< MATRIX >::disable_resolve(), oomph::HelmholtzGMRESMG< MATRIX >::solve(), oomph::HelmholtzFGMRESMG< MATRIX >::solve(), oomph::HelmholtzFGMRESMG< MATRIX >::~HelmholtzFGMRESMG(), and oomph::HelmholtzGMRESMG< MATRIX >::~HelmholtzGMRESMG().

◆ complex_matrix_multiplication()

template<typename MATRIX >

void oomph::HelmholtzGMRESMG< MATRIX >::complex_matrix_multiplication	(	Vector< CRDoubleMatrix * > const	matrices_pt,
		const Vector< DoubleVector > &	x,
		Vector< DoubleVector > &	soln
	)

inlineprotected

Helper function to calculate a complex matrix-vector product. Assumes the matrix has been provided as a Vector of length two; the first entry containing the real part of the system matrix and the second entry containing the imaginary part

     {
 #ifdef PARANOID
       // PARANOID check - Make sure the input matrix has the right size
       if (matrices_pt.size() != 2)
       {
         // Create an output stream
         std::ostringstream error_message_stream;
  
         // Create the error message
         error_message_stream << "Can only deal with two matrices. You have "
                              << matrices_pt.size() << " matrices." << std::endl;
  
         // Throw an error
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
       // PARANOID check - Make sure the vector x has the right size
       if (x.size() != 2)
       {
         // Create an output stream
         std::ostringstream error_message_stream;
  
         // Create the error message
         error_message_stream
           << "Can only deal with two input vectors. You have " << x.size()
           << " vectors." << std::endl;
  
         // Throw an error
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
       // PARANOID check - Make sure the vector soln has the right size
       if (soln.size() != 2)
       {
         // Create an output stream
         std::ostringstream error_message_stream;
  
         // Create the error message
         error_message_stream
           << "Can only deal with two output vectors. You have " << soln.size()
           << " output vectors." << std::endl;
  
         // Throw an error
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       // NOTE: We assume all vectors have been distributed at this point but
       // code can be written at a later time to build the vectors if they're
       // not already built.
  
       //-----------------------------------------------------------------------
       // Suppose we have a complex matrix, A, and two complex vectors, x and
       // soln. We wish to compute the product A*x=soln (note, * does not mean
       // we are using complex conjugates here, it is simply used to indicate
       // a multiplication). To do this we must make use of the fact that we
       // possess the real and imaginary separately. As a result, it is computed
       // using:
       //           soln = A*x,
       //                = (A_r + i*A_c)*(x_r + i*x_c),
       //                = [A_r*x_r - A_c*x_c] + i*[A_r*x_c + A_c*x_r],
       // ==> real(soln) = A_r*x_r - A_c*x_c,
       //   & imag(soln) = A_r*x_c + A_c*x_r,
       // where the subscripts _r and _c are used to identify the real and
       // imaginary part, respectively.
       //-----------------------------------------------------------------------
  
       // Store the value of A_r*x_r in the real part of soln
       matrices_pt[0]->multiply(x[0], soln[0]);
  
       // Store the value of A_r*x_c in the imaginary part of soln
       matrices_pt[0]->multiply(x[1], soln[1]);
  
       // Create a temporary vector
       DoubleVector temp(Matrices_storage_pt[0]->distribution_pt(), 0.0);
  
       // Calculate the value of A_c*x_c
       matrices_pt[1]->multiply(x[1], temp);
  
       // Subtract the value of temp from soln[0] to get the real part of soln
       soln[0] -= temp;
  
       // Calculate the value of A_c*x_r
       matrices_pt[1]->multiply(x[0], temp);
  
       // Add the value of temp to soln[1] to get the imaginary part of soln
       soln[1] += temp;
     } // End of complex_matrix_multiplication

References oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::HelmholtzGMRESMG< MATRIX >::Matrices_storage_pt, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and plotDoE::x.

◆ disable_resolve()

template<typename MATRIX >

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

inlinevirtual

Overload disable resolve so that it cleans up memory too.

Reimplemented from oomph::LinearSolver.

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

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

◆ generate_plane_rotation()

template<typename MATRIX >

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

inlineprotected

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

\[ \begin{bmatrix} \overline{\cos\theta} & \overline{\sin\theta} \\ -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \\ dy \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix}, \]

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

\[ \cos\theta&=\dfrac{dx}{\sqrt{|dx|^2+|dy|^2}}, \]

and

\[ \sin\theta&=\dfrac{dy}{\sqrt{|dx|^2+|dy|^2}}. \]

Taken from: Saad Y."Iterative methods for sparse linear systems", p.193. We also check to see that sn is always a real (nonnegative) number. See pp.193-194 for an explanation.

     {
       // 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 original 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)) 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 (std::abs(dy) > std::abs(dx))
       {
         // Since |dy|>|dx| calculate the ratio |dx|/|dy|
         std::complex<double> temp = dx / dy;
  
         // Calculate the value of sin(theta) using:
         //          sin(theta)=dy/sqrt(pow(|dx|,2)+pow(|dy|,2))
         //                    =(dy/|dy|)/sqrt(1+pow(|dx|/|dy|,2)).
         sn = (dy / std::abs(dy)) / sqrt(1.0 + pow(std::abs(temp), 2.0));
  
         // Calculate the value of cos(theta) using:
         //            cos(theta)=dx/sqrt(pow(|dy|,2)+pow(|dx|,2))
         //                      =(dx/|dy|)/sqrt(1+pow(|dx|/|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
         std::complex<double> temp = dy / dx;
  
         // Calculate the value of cos(theta) using:
         //          cos(theta)=dx/sqrt(pow(|dx|,2)+pow(|dy|,2))
         //                    =(dx/|dx|)/sqrt(1+pow(|dy|/|dx|,2)).
         cs = (dx / std::abs(dx)) / sqrt(1.0 + pow(std::abs(temp), 2.0));
  
         // Calculate the value of sin(theta) using:
         //            sin(theta)=dy/sqrt(pow(|dx|,2)+pow(|dy|,2))
         //                      =(dy/|dx|)/sqrt(1+pow(|dy|/|dx|,2))
         //                      =(dy/dx)*cos(theta).
         sn = temp * cs;
       }
  
       // Set the tolerance for sin(theta)
       double tolerance = 1.0e-15;
  
       // Make sure sn is real and nonnegative (it should be!)
       if ((std::fabs(sn.imag()) > tolerance) || (sn.real() < 0))
       {
         // Create an output stream
         std::ostringstream error_message_stream;
  
         // Create an error message
         error_message_stream << "The value of sin(theta) is not real "
                              << "and/or nonnegative. Value is: " << sn
                              << std::endl;
  
         // Throw an error
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     } // End of generate_plane_rotation

References abs(), boost::multiprecision::fabs(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Eigen::bfloat16_impl::pow(), sqrt(), and oomph::IterativeLinearSolver::tolerance().

◆ iterations()

template<typename MATRIX >

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

inlinevirtual

Number of iterations taken.

Implements oomph::IterativeLinearSolver.

     {
       return Iterations;
     }

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

◆ operator=()

template<typename MATRIX >

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

delete

Broken assignment operator.

◆ preconditioner_solve()

template<typename MATRIX >

void oomph::HelmholtzGMRESMG< MATRIX >::preconditioner_solve	(	const DoubleVector &	r,
		DoubleVector &	z
	)

inlinevirtual

Implementation of the pure virtual base class function. The function has been broken because this is meant to be used as a linear solver

Implements oomph::Preconditioner.

     {
       // Open an output stream
       std::ostringstream error_message_stream;
  
       // Create an error message
       error_message_stream << "Preconditioner_solve(...) is broken. "
                            << "HelmholtzGMRESMG is only meant to be used as "
                            << "a linear solver.\n";
  
       // Throw the error message
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ resolve()

template<typename MATRIX >

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

inlinevirtual

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 ((Matrices_storage_pt[0] == 0) || (Matrices_storage_pt[1] == 0))
       {
         throw OomphLibError("No matrix was stored -- cannot re-solve",
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       // Set up a dummy matrix. As we're resolving this won't be used in
       // solve_helper but we need to pass a matrix in to fill the input.
       // The matrices used in the calculations have already been stored
       CRDoubleMatrix* matrix_pt = new CRDoubleMatrix;
  
       // Call the helper function
       solve_helper(matrix_pt, rhs, result);
  
       // Delete the matrix
       delete matrix_pt;
  
       // Make it a null pointer
       matrix_pt = 0;
  
       // Reset re-solving flag
       Resolving = false;
     }

References oomph::HelmholtzGMRESMG< MATRIX >::Matrices_storage_pt, oomph::BlockPreconditioner< MATRIX >::matrix_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::HelmholtzGMRESMG< MATRIX >::Resolving, and oomph::HelmholtzGMRESMG< MATRIX >::solve_helper().

◆ set_preconditioner_LHS()

template<typename MATRIX >

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

inline

Set left preconditioning (the default)

     {
       Preconditioner_LHS = true;
     }

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

◆ set_preconditioner_RHS()

template<typename MATRIX >

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

inline

Enable right preconditioning.

     {
       Preconditioner_LHS = false;
     }

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

Referenced by PMLHelmholtzMGProblem< ELEMENT >::set_gmres_multigrid_solver().

◆ setup()

template<typename MATRIX >

void oomph::HelmholtzGMRESMG< MATRIX >::setup ( )

inlinevirtual

Implementation of the pure virtual base class function. This accompanies the preconditioner_solve function and so is also broken

Implements oomph::Preconditioner.

     {
       // Open an output stream
       std::ostringstream error_message_stream;
  
       // Create an error message
       error_message_stream << "This function is broken. HelmholtzGMRESMG is "
                            << "only meant to be used as a linear solver.\n";
  
       // Throw the error message
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ solve() [1/3]

template<typename MATRIX >

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

inlinevirtual

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

Reimplemented from oomph::LinearSolver.

     {
       // Open an output stream
       std::ostringstream error_message_stream;
  
       // Create an error message
       error_message_stream
         << "This function is broken. The block preconditioner "
         << "needs access to the underlying mesh.\n";
  
       // Throw the error message
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ solve() [2/3]

template<typename MATRIX >

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

inlinevirtual

Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the linear system Call the broken base-class version. If you want this, please implement it

Reimplemented from oomph::LinearSolver.

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

References oomph::BlockPreconditioner< MATRIX >::matrix_pt(), and oomph::LinearSolver::solve().

◆ solve() [3/3]

template<typename MATRIX >

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

inlinevirtual

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.

Reimplemented in oomph::HelmholtzFGMRESMG< MATRIX >.

     {
 #ifdef OOMPH_HAS_MPI
       // Make sure that this is running in serial. Can't guarantee it'll
       // work when the problem is distributed over several processors
       if (MPI_Helpers::communicator_pt()->nproc() > 1)
       {
         // Throw a warning
         OomphLibWarning("Can't guarantee the MG solver will work in parallel!",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       // Find # of degrees of freedom (variables)
       unsigned n_dof = problem_pt->ndof();
  
       // Initialise timer
       double t_start = TimingHelpers::timer();
  
       // We're not re-solving
       Resolving = false;
  
       // Get rid of any previously stored data
       clean_up_memory();
  
       // Grab the communicator from the MGProblem object and assign it
       this->set_comm_pt(problem_pt->communicator_pt());
  
       // Setup the distribution
       LinearAlgebraDistribution dist(
         problem_pt->communicator_pt(), n_dof, false);
  
       // Build the internal distribution in this way because both the
       // IterativeLinearSolver and BlockPreconditioner class have base-
       // class subobjects of type oomph::DistributableLinearAlgebraObject
       IterativeLinearSolver::build_distribution(dist);
  
       // Get Jacobian matrix in format specified by template parameter
       // and nonlinear residual vector
       MATRIX* matrix_pt = new MATRIX;
       DoubleVector f;
       if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
       {
         if (dynamic_cast<CRDoubleMatrix*>(matrix_pt) != 0)
         {
           dynamic_cast<CRDoubleMatrix*>(matrix_pt)->build(
             IterativeLinearSolver::distribution_pt());
           f.build(IterativeLinearSolver::distribution_pt(), 0.0);
         }
       }
  
       // Get the Jacobian and residuals vector
       problem_pt->get_jacobian(f, *matrix_pt);
  
       // We've made the matrix, we can delete it...
       Matrix_can_be_deleted = true;
  
       // Replace the current matrix used in Preconditioner by the new matrix
       this->set_matrix_pt(matrix_pt);
  
       // The preconditioner works with one mesh; set it! Since we only use
       // the block preconditioner on the finest level, we use the mesh from
       // that level
       this->set_nmesh(1);
  
       // Elements in actual pml layer are trivially wrapped versions of
       // their bulk counterparts. Technically they are different elements
       // so we have to allow different element types
       bool allow_different_element_types_in_mesh = true;
       this->set_mesh(
         0, problem_pt->mesh_pt(), allow_different_element_types_in_mesh);
  
       // Set up the generic block look up scheme
       this->block_setup();
  
       // Extract the number of blocks.
       unsigned nblock_types = this->nblock_types();
  
 #ifdef PARANOID
       // PARANOID check - there must only be two block types
       if (nblock_types != 2)
       {
         // Create the error message
         std::stringstream tmp;
         tmp << "There are supposed to be two block types.\nYours has "
             << nblock_types << std::endl;
  
         // Throw an error
         throw OomphLibError(
           tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       // Resize the storage for the system matrices
       Matrices_storage_pt.resize(2, 0);
  
       // Loop over the rows of the block matrix
       for (unsigned i = 0; i < nblock_types; i++)
       {
         // Fix the column index
         unsigned j = 0;
  
         // Create new CRDoubleMatrix objects
         Matrices_storage_pt[i] = new CRDoubleMatrix;
  
         // Extract the required blocks, i.e. the first column
         this->get_block(i, j, *Matrices_storage_pt[i]);
       }
  
       // Doc time for setup
       double t_end = TimingHelpers::timer();
       Jacobian_setup_time = t_end - t_start;
  
       if (Doc_time)
       {
         oomph_info << "\nTime for setup of block Jacobian [sec]: "
                    << Jacobian_setup_time << std::endl;
       }
  
       // Call linear algebra-style solver
       // If the result distribution is wrong, then redistribute
       // before the solve and return to original distribution
       // afterwards
       if ((!(*result.distribution_pt() ==
              *IterativeLinearSolver::distribution_pt())) &&
           (result.built()))
       {
         // Make a distribution object
         LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
  
         // Build the result vector distribution
         result.build(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Solve the problem
         solve_helper(matrix_pt, f, result);
  
         // Redistribute the vector
         result.redistribute(&temp_global_dist);
       }
       // Otherwise just solve
       else
       {
         // Solve
         solve_helper(matrix_pt, f, result);
       }
  
       // Kill matrix unless it's still required for resolve
       if (!Enable_resolve)
       {
         // Clean up anything left in memory
         clean_up_memory();
       }
     } // End of solve

References oomph::BlockPreconditioner< MATRIX >::block_setup(), oomph::DoubleVector::build(), oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::built(), oomph::HelmholtzGMRESMG< MATRIX >::clean_up_memory(), oomph::MPI_Helpers::communicator_pt(), oomph::Problem::communicator_pt(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Doc_time, oomph::LinearSolver::Enable_resolve, f(), oomph::BlockPreconditioner< MATRIX >::get_block(), oomph::Problem::get_jacobian(), i, j, oomph::IterativeLinearSolver::Jacobian_setup_time, oomph::HelmholtzGMRESMG< MATRIX >::Matrices_storage_pt, oomph::HelmholtzGMRESMG< MATRIX >::Matrix_can_be_deleted, oomph::BlockPreconditioner< MATRIX >::matrix_pt(), oomph::Problem::mesh_pt(), oomph::BlockPreconditioner< MATRIX >::nblock_types(), oomph::Problem::ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::DoubleVector::redistribute(), oomph::HelmholtzGMRESMG< MATRIX >::Resolving, oomph::Preconditioner::set_comm_pt(), oomph::Preconditioner::set_matrix_pt(), oomph::BlockPreconditioner< MATRIX >::set_mesh(), oomph::BlockPreconditioner< MATRIX >::set_nmesh(), oomph::HelmholtzGMRESMG< MATRIX >::solve_helper(), oomph::TimingHelpers::timer(), and tmp.

◆ solve_helper()

template<typename MATRIX >

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

protected

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

   {
     // Set the number of dof types (real and imaginary for this solver)
     unsigned n_dof_types = this->ndof_types();
  
 #ifdef PARANOID
     // This only works for 2 dof types
     if (n_dof_types != 2)
     {
       // Create an output stream
       std::stringstream error_message_stream;
  
       // Create the error message
       error_message_stream << "This preconditioner only works for problems "
                            << "with 2 dof types\nYours has " << n_dof_types;
  
       // Throw the error message
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Get the number of dofs (note, the total number of dofs in the problem
     // is 2*n_row but if the constituent vectors and matrices were stored in
     // complex objects there would only be (n_row) rows so we use that)
     unsigned n_row = Matrices_storage_pt[0]->nrow();
  
     // Make sure Max_iter isn't greater than n_dof. The user cannot use this
     // many iterations when using Krylov subspace methods
     if (Max_iter > n_row)
     {
       // Create an output string stream
       std::ostringstream error_message_stream;
  
       // Create the error message
       error_message_stream << "The maximum number of iterations cannot exceed "
                            << "the number of rows in the problem."
                            << "\nMaximum number of iterations: " << Max_iter
                            << "\nNumber of rows: " << n_row << std::endl;
  
       // Throw the error message
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
 #ifdef PARANOID
     // Loop over the real and imaginary parts
     for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
     {
       // PARANOID check that if the matrix is distributable then it should not
       // be then it should not be distributed
       if (dynamic_cast<DistributableLinearAlgebraObject*>(
             Matrices_storage_pt[dof_type]) != 0)
       {
         if (dynamic_cast<DistributableLinearAlgebraObject*>(
               Matrices_storage_pt[dof_type])
               ->distributed())
         {
           std::ostringstream error_message_stream;
           error_message_stream << "The matrix must not be distributed.";
           throw OomphLibError(error_message_stream.str(),
                               OOMPH_CURRENT_FUNCTION,
                               OOMPH_EXCEPTION_LOCATION);
         }
       }
     }
     // PARANOID check that this rhs distribution is setup
     if (!rhs.built())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The rhs vector distribution must be setup.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // PARANOID check that the rhs has the right number of global rows
     if (rhs.nrow() != 2 * n_row)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "RHS does not have the same dimension as the"
                            << " linear system";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // PARANOID check that the rhs is not distributed
     if (rhs.distribution_pt()->distributed())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The rhs vector must not be distributed.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // PARANOID check that if the result is setup it matches the distribution
     // of the rhs
     if (solution.built())
     {
       if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
       {
         std::ostringstream error_message_stream;
         error_message_stream << "If the result distribution is setup then it "
                              << "must be the same as the rhs distribution";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     } // if (solution[dof_type].built())
 #endif
  
     // Set up the solution distribution if it's not already distributed
     if (!solution.built())
     {
       // Build the distribution
       solution.build(IterativeLinearSolver::distribution_pt(), 0.0);
     }
     // Otherwise initialise all entries to zero
     else
     {
       // Initialise the entries in the k-th vector in solution to zero
       solution.initialise(0.0);
     }
  
     // Create a vector of DoubleVectors (this is a distributed vector so we have
     // to create two separate DoubleVector objects to cope with the arithmetic)
     Vector<DoubleVector> block_solution(n_dof_types);
  
     // Create a vector of DoubleVectors (this is a distributed vector so we have
     // to create two separate DoubleVector objects to cope with the arithmetic)
     Vector<DoubleVector> block_rhs(n_dof_types);
  
     // Build the distribution of both vectors
     for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
     {
       // Build the distribution of the k-th constituent vector
       block_solution[dof_type].build(this->block_distribution_pt(dof_type),
                                      0.0);
  
       // Build the distribution of the k-th constituent vector
       block_rhs[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
     }
  
     // Grab the solution vector in block form
     this->get_block_vectors(solution, block_solution);
  
     // Grab the RHS vector in block form
     this->get_block_vectors(rhs, block_rhs);
  
     // Start the solver timer
     double t_start = TimingHelpers::timer();
  
     // Storage for the relative residual
     double resid;
  
     // Initialise vectors (since these are not distributed vectors we template
     // one vector by the type std::complex<double> instead of using two vectors,
     // each templated by the type double
  
     // Vector, s, used in the minimisation problem: J(y)=min||s-R_m*y||
     // [see Saad Y."Iterative methods for sparse linear systems", p.176.]
     Vector<std::complex<double>> s(n_row + 1, std::complex<double>(0.0, 0.0));
  
     // Vector to store the value of cos(theta) when using the Givens rotation
     Vector<std::complex<double>> cs(n_row + 1, std::complex<double>(0.0, 0.0));
  
     // Vector to store the value of sin(theta) when using the Givens rotation
     Vector<std::complex<double>> sn(n_row + 1, std::complex<double>(0.0, 0.0));
  
     // Create a vector of DoubleVectors (this is a distributed vector so we have
     // to create two separate DoubleVector objects to cope with the arithmetic)
     Vector<DoubleVector> block_w(n_dof_types);
  
     // Build the distribution of both vectors
     for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
     {
       // Build the distribution of the k-th constituent vector
       block_w[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
     }
  
     // Set up the preconditioner only if we're not re-solving
     if (!Resolving)
     {
       // Only set up the preconditioner before solve if required
       if (Setup_preconditioner_before_solve)
       {
         // Set up the preconditioner from the Jacobian matrix
         double t_start_prec = TimingHelpers::timer();
  
         // Use the setup function in the Preconditioner class
         preconditioner_pt()->setup(dynamic_cast<MATRIX*>(matrix_pt));
  
         // Doc time for setup of preconditioner
         double t_end_prec = TimingHelpers::timer();
         Preconditioner_setup_time = t_end_prec - t_start_prec;
  
         // If time documentation is enabled
         if (Doc_time)
         {
           // Output the time taken
           oomph_info << "Time for setup of preconditioner [sec]: "
                      << Preconditioner_setup_time << std::endl;
         }
       }
     }
     else
     {
       // If time documentation is enabled
       if (Doc_time)
       {
         // Notify the user
         oomph_info << "Setup of preconditioner is bypassed in resolve mode"
                    << std::endl;
       }
     } // if (!Resolving) else
  
     // Create a vector of DoubleVectors to store the RHS of b-Jx=Mr. We assume
     // both x=0 and that a preconditioner is not applied by which we deduce b=r
     Vector<DoubleVector> block_r(n_dof_types);
  
     // Build the distribution of both vectors
     for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
     {
       // Build the distribution of the k-th constituent vector
       block_r[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
     }
  
     // If we're using LHS preconditioning solve b-Jx=Mr for r (assumes x=0)
     // so calculate r=M^{-1}b otherwise set r=b (RHS prec.)
     if (Preconditioner_LHS)
     {
       // Create a vector of the same size as rhs
       DoubleVector r(IterativeLinearSolver::distribution_pt(), 0.0);
  
       // Copy the vectors in r to full_r
       this->return_block_vectors(block_r, r);
  
       // Use the preconditioner
       preconditioner_pt()->preconditioner_solve(rhs, r);
  
       // Copy the vector full_r into the vectors in r
       this->get_block_vectors(r, block_r);
     }
     else
     {
       // Store the value of b (the RHS vector) in r
       for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
       {
         // Copy the entries of rhs into r
         block_r[dof_type] = block_rhs[dof_type];
       }
     } // if(Preconditioner_LHS)
  
     // Calculate the norm of the real part of r
     double norm_r = block_r[0].norm();
  
     // Calculate the norm of the imaginary part of r
     double norm_c = block_r[1].norm();
  
     // Compute norm(r)
     double normb = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
  
     // Set the value of beta (the initial residual)
     double beta = normb;
  
     // Compute the initial relative residual. If the entries of the RHS vector
     // are all zero then set normb equal to one. This is because we divide the
     // value of the norm at later stages by normb and dividing by zero is not
     // definied
     if (normb == 0.0)
     {
       // Set the value of normb
       normb = 1.0;
     }
  
     // Calculate the ratio between the initial norm and the current norm.
     // Since we haven't completed an iteration yet this will simply be one
     // unless normb was zero, in which case resid will have value zero
     resid = beta / normb;
  
     // If required, will document convergence history to screen or file (if
     // stream open)
     if (Doc_convergence_history)
     {
       // If an output file which is open isn't provided then output to screen
       if (!Output_file_stream.is_open())
       {
         // Output the residual value to the screen
         oomph_info << 0 << " " << resid << std::endl;
       }
       // Otherwise, output to file
       else
       {
         // Output the residual value to file
         Output_file_stream << 0 << " " << resid << std::endl;
       }
     } // if (Doc_convergence_history)
  
     // If the GMRES algorithm converges immediately
     if (resid <= Tolerance)
     {
       // If time documentation is enabled
       if (Doc_time)
       {
         // Notify the user
         oomph_info << "GMRES converged immediately. Normalised residual norm: "
                    << resid << std::endl;
       }
  
       // Finish running the solver
       return;
     } // if (resid<=Tolerance)
  
     // Initialise a vector of orthogonal basis vectors
     Vector<Vector<DoubleVector>> block_v;
  
     // Resize the number of vectors needed
     block_v.resize(n_row + 1);
  
     // Resize each Vector of DoubleVectors to store the real and imaginary
     // part of a given vector
     for (unsigned dof_type = 0; dof_type < n_row + 1; dof_type++)
     {
       // Create two DoubleVector objects in each Vector
       block_v[dof_type].resize(n_dof_types);
     }
  
     // Initialise the upper hessenberg matrix. Since we are not using
     // distributed vectors here, the algebra is best done using entries
     // of the type std::complex<double>. NOTE: For implementation purposes
     // the upper Hessenberg matrix indices are swapped so the matrix is
     // effectively transposed
     Vector<Vector<std::complex<double>>> hessenberg(n_row + 1);
  
     // Build the zeroth basis vector
     for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
     {
       // Build the k-th part of the zeroth vector. Here k=0 and k=1 correspond
       // to the real and imaginary part of the zeroth vector, respectively
       block_v[0][dof_type].build(this->block_distribution_pt(dof_type), 0.0);
     }
  
     // Loop over the real and imaginary parts of v
     for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
     {
       // For fast access
       double* v0_pt = block_v[0][dof_type].values_pt();
  
       // For fast access
       const double* block_r_pt = block_r[dof_type].values_pt();
  
       // Set the zeroth basis vector v[0] to r/beta
       for (unsigned i = 0; i < n_row; i++)
       {
         // Assign the i-th entry of the zeroth basis vector
         v0_pt[i] = block_r_pt[i] / beta;
       }
     } // for (unsigned k=0;k<n_dof_types;k++)
  
     // Set the first entry in the minimisation problem RHS vector (is meant
     // to the vector beta*e_1 initially, where e_1 is the unit vector with
     // one in its first entry)
     s[0] = beta;
  
     // Compute the next step of the iterative scheme
     for (unsigned j = 0; j < Max_iter; j++)
     {
       // Resize the next column of the upper hessenberg matrix
       hessenberg[j].resize(j + 2, std::complex<double>(0.0, 0.0));
  
       // Calculate w=M^{-1}(Jv[j]) (LHS prec.) or w=JM^{-1}v (RHS prec.)
       {
         // Create a temporary vector
         DoubleVector vj(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Create a temporary vector
         DoubleVector temp(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Create a temporary vector
         DoubleVector w(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Create a temporary vector of DoubleVectors
         Vector<DoubleVector> block_temp(2);
  
         // Create two DoubleVectors
         for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
         {
           block_temp[dof_type].build(this->block_distribution_pt(dof_type),
                                      0.0);
         }
  
         // If we're using LHS preconditioning
         if (Preconditioner_LHS)
         {
           // Solve Jv[j]=Mw for w. Note, we cannot use inbuilt complex matrix
           // algebra here as we're using distributed vectors
           complex_matrix_multiplication(
             Matrices_storage_pt, block_v[j], block_temp);
  
           // Copy block_temp into temp
           this->return_block_vectors(block_temp, temp);
  
           // Copy block_w into w
           this->return_block_vectors(block_w, w);
  
           // Apply the preconditioner
           this->preconditioner_pt()->preconditioner_solve(temp, w);
  
           // Copy w into block_w
           this->get_block_vectors(w, block_w);
         }
         // If we're using RHS preconditioning
         else
         {
           // Copy the real and imaginary part of v[j] into one vector, vj
           this->return_block_vectors(block_v[j], vj);
  
           // Use w=JM^{-1}v by saad p270
           this->preconditioner_pt()->preconditioner_solve(vj, temp);
  
           // Copy w into block_w
           this->get_block_vectors(temp, block_temp);
  
           // Solve Jv[j] = Mw for w. Note, we cannot use inbuilt complex matrix
           // algebra here as we're using distributed vectors
           complex_matrix_multiplication(
             Matrices_storage_pt, block_temp, block_w);
         }
       } // Calculate w=M^{-1}(Jv[j]) (LHS prec.) or w=JM^{-1}v (RHS prec.)
  
       // For fast access
       double* block_w_r_pt = block_w[0].values_pt();
  
       // For fast access
       double* block_w_c_pt = block_w[1].values_pt();
  
       // Loop over all of the entries on and above the principal subdiagonal of
       // the Hessenberg matrix in the j-th column (remembering that
       // the indices of the upper Hessenberg matrix are swapped for the purpose
       // of implementation)
       for (unsigned i = 0; i < j + 1; i++)
       {
         // For fast access
         const double* vi_r_pt = block_v[i][0].values_pt();
  
         // For fast access
         const double* vi_c_pt = block_v[i][1].values_pt();
  
         // Loop over the entries of v and w
         for (unsigned k = 0; k < n_row; k++)
         {
           // Store the appropriate entry in v as a complex value
           std::complex<double> complex_v(vi_r_pt[k], vi_c_pt[k]);
  
           // Store the appropriate entry in w as a complex value
           std::complex<double> complex_w(block_w_r_pt[k], block_w_c_pt[k]);
  
           // Update the value of H(i,j) noting we're computing a complex
           // inner product here (the ordering is very important here!)
           hessenberg[j][i] += std::conj(complex_v) * complex_w;
         }
  
         // Orthonormalise w against all previous orthogonal vectors, v_i by
         // looping over its entries and updating them
         for (unsigned k = 0; k < n_row; k++)
         {
           // Update the real part of the k-th entry of w
           block_w_r_pt[k] -= (hessenberg[j][i].real() * vi_r_pt[k] -
                               hessenberg[j][i].imag() * vi_c_pt[k]);
  
           // Update the imaginary part of the k-th entry of w
           block_w_c_pt[k] -= (hessenberg[j][i].real() * vi_c_pt[k] +
                               hessenberg[j][i].imag() * vi_r_pt[k]);
         }
       } // for (unsigned i=0;i<j+1;i++)
  
       // Calculate the norm of the real part of w
       norm_r = block_w[0].norm();
  
       // Calculate the norm of the imaginary part of w
       norm_c = block_w[1].norm();
  
       // Calculate the norm of the vector w using norm_r and norm_c and assign
       // its value to the appropriate entry in the Hessenberg matrix
       hessenberg[j][j + 1] = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
  
       // Build the (j+1)-th basis vector
       for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
       {
         // Build the k-th part of the zeroth vector. Here k=0 and k=1 correspond
         // to the real and imaginary part of the zeroth vector, respectively
         block_v[j + 1][dof_type].build(this->block_distribution_pt(dof_type),
                                        0.0);
       }
  
       // Check if the value of hessenberg[j][j+1] is zero. If it
       // isn't then we update the next entries in v
       if (hessenberg[j][j + 1] != 0.0)
       {
         // For fast access
         double* v_r_pt = block_v[j + 1][0].values_pt();
  
         // For fast access
         double* v_c_pt = block_v[j + 1][1].values_pt();
  
         // For fast access
         const double* block_w_r_pt = block_w[0].values_pt();
  
         // For fast access
         const double* block_w_c_pt = block_w[1].values_pt();
  
         // Notice, the value of H(j,j+1), as calculated above, is clearly a real
         // number. As such, calculating the division
         //                      v_{j+1}=w_{j}/h_{j+1,j},
         // here is simple, i.e. we don't need to worry about cross terms in the
         // algebra. To avoid computing h_{j+1,j} several times we precompute it
         double h_subdiag_val = hessenberg[j][j + 1].real();
  
         // Loop over the entries of the new orthogonal vector and set its values
         for (unsigned k = 0; k < n_row; k++)
         {
           // The i-th entry of the real component is given by
           v_r_pt[k] = block_w_r_pt[k] / h_subdiag_val;
  
           // Similarly, the i-th entry of the imaginary component is given by
           v_c_pt[k] = block_w_c_pt[k] / h_subdiag_val;
         }
       }
       // Otherwise, we have to jump to the next part of the algorithm; if
       // the value of hessenberg[j][j+1] is zero then the norm of the latest
       // orthogonal vector is zero. This is only possible if the entries
       // in w are all zero. As a result, the Krylov space of A and r_0 has
       // been spanned by the previously calculated orthogonal vectors
       else
       {
         // Book says "Set m=j and jump to step 11" (p.172)...
         // Do something here!
         oomph_info << "Subdiagonal Hessenberg entry is zero. "
                    << "Do something here..." << std::endl;
       } // if (hessenberg[j][j+1]!=0.0)
  
       // Loop over the entries in the Hessenberg matrix and calculate the
       // entries of the Givens rotation matrices
       for (unsigned k = 0; k < j; k++)
       {
         // Apply the plane rotation to all of the previous entries in the
         // (j)-th column (remembering the indexing is reversed)
         apply_plane_rotation(
           hessenberg[j][k], hessenberg[j][k + 1], cs[k], sn[k]);
       }
  
       // Now calculate the entries of the latest Givens rotation matrix
       generate_plane_rotation(
         hessenberg[j][j], hessenberg[j][j + 1], cs[j], sn[j]);
  
       // Apply the plane rotation using the newly calculated entries
       apply_plane_rotation(
         hessenberg[j][j], hessenberg[j][j + 1], cs[j], sn[j]);
  
       // Apply a plane rotation to the corresponding entry in the vector
       // s used in the minimisation problem, J(y)=min||s-R_m*y||
       apply_plane_rotation(s[j], s[j + 1], cs[j], sn[j]);
  
       // Compute current residual using equation (6.42) in Saad Y, "Iterative
       // methods for sparse linear systems", p.177.]. Note, since s has complex
       // entries we have to use std::abs instead of std::fabs
       beta = std::abs(s[j + 1]);
  
       // 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 an output file which is open isn't provided then output to screen
         if (!Output_file_stream.is_open())
         {
           // Output the residual value to the screen
           oomph_info << j + 1 << " " << resid << std::endl;
         }
         // Otherwise, output to file
         else
         {
           // Output the residual value to file
           Output_file_stream << j + 1 << " " << resid << std::endl;
         }
       } // if (Doc_convergence_history)
  
       // If the required tolerance has been met
       if (resid < Tolerance)
       {
         // Store the number of iterations taken
         Iterations = j + 1;
  
         // Update the result vector using the result, x=x_0+V_m*y (where V_m
         // is given by v here)
         update(j, hessenberg, s, block_v, block_solution);
  
         // Copy the vectors in block_solution to solution
         this->return_block_vectors(block_solution, solution);
  
         // If time documentation was enabled
         if (Doc_time)
         {
           // Output the current normalised residual norm
           oomph_info << "\nGMRES converged (1). Normalised residual norm: "
                      << resid << std::endl;
  
           // Output the number of iterations it took for convergence
           oomph_info << "Number of iterations to convergence: " << j + 1 << "\n"
                      << std::endl;
         }
  
         // Stop the timer
         double t_end = TimingHelpers::timer();
  
         // Calculate the time taken to calculate the solution
         Solution_time = t_end - t_start;
  
         // If time documentation was enabled
         if (Doc_time)
         {
           // Output the time taken to solve the problem using GMRES
           oomph_info << "Time for solve with GMRES [sec]: " << Solution_time
                      << std::endl;
         }
  
         // As we've met the tolerance for the solver and everything that should
         // be documented, has been, finish using the solver
         return;
       } // if (resid<Tolerance)
     } // for (unsigned j=0;j<Max_iter;j++)
  
     // Store the number of iterations taken
     Iterations = Max_iter;
  
     // Only update if we actually did something
     if (Max_iter > 0)
     {
       // Update the result vector using the result, x=x_0+V_m*y (where V_m
       // is given by v here)
       update(Max_iter - 1, hessenberg, s, block_v, block_solution);
  
       // Copy the vectors in block_solution to solution
       this->return_block_vectors(block_solution, solution);
     }
  
     // Solve Mr=(b-Jx) for r
     {
       // Create a temporary vector of DoubleVectors
       Vector<DoubleVector> block_temp(2);
  
       // Create two DoubleVectors
       for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
       {
         // Build the distribution of the (dof_type)-th vector
         block_temp[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
       }
  
       // Calculate the value of Jx
       complex_matrix_multiplication(
         Matrices_storage_pt, block_solution, block_temp);
  
       // Get the values pointer of the vector (real)
       double* block_temp_r_pt = block_temp[0].values_pt();
  
       // Get the values pointer of the vector (imaginary)
       double* block_temp_c_pt = block_temp[1].values_pt();
  
       // Get the values pointer of the RHS vector (real)
       const double* block_rhs_r_pt = block_rhs[0].values_pt();
  
       // Get the values pointer of the RHS vector (imaginary)
       const double* block_rhs_c_pt = block_rhs[1].values_pt();
  
       // Loop over the dofs
       for (unsigned i = 0; i < n_row; i++)
       {
         // Calculate b-Jx (real)
         block_temp_r_pt[i] = block_rhs_r_pt[i] - block_temp_r_pt[i];
  
         // Calculate b-Jx (imaginary)
         block_temp_c_pt[i] = block_rhs_c_pt[i] - block_temp_c_pt[i];
       }
  
       // If we're using LHS preconditioning
       if (Preconditioner_LHS)
       {
         // Create a temporary DoubleVectors
         DoubleVector temp(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Create a vector of the same size as rhs
         DoubleVector r(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Copy the vectors in r to full_r
         this->return_block_vectors(block_temp, temp);
  
         // Copy the vectors in r to full_r
         this->return_block_vectors(block_r, r);
  
         // Apply the preconditioner
         preconditioner_pt()->preconditioner_solve(temp, r);
       }
     }
  
     // Compute the current residual
     beta = 0.0;
  
     // Get access to the values pointer (real)
     norm_r = block_r[0].norm();
  
     // Get access to the values pointer (imaginary)
     norm_c = block_r[1].norm();
  
     // Calculate the full norm
     beta = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
  
     // Calculate the relative residual
     resid = beta / normb;
  
     // If the relative residual lies within tolerance
     if (resid < Tolerance)
     {
       // If time documentation is enabled
       if (Doc_time)
       {
         // Notify the user
         oomph_info << "\nGMRES converged (2). Normalised residual norm: "
                    << resid
                    << "\nNumber of iterations to convergence: " << Iterations
                    << "\n"
                    << std::endl;
       }
  
       // End the timer
       double t_end = TimingHelpers::timer();
  
       // Calculate the time taken for the solver
       Solution_time = t_end - t_start;
  
       // If time documentation is enabled
       if (Doc_time)
       {
         oomph_info << "Time for solve with GMRES [sec]: " << Solution_time
                    << std::endl;
       }
       return;
     }
  
     // Otherwise GMRES failed convergence
     oomph_info << "\nGMRES did not converge to required tolerance! "
                << "\nReturning with normalised residual norm: " << resid
                << "\nafter " << Max_iter << " iterations.\n"
                << std::endl;
  
     // Throw an error if requested
     if (Throw_error_after_max_iter)
     {
       std::string err = "Solver failed to converge and you requested an error";
       err += " on convergence failures.";
       throw OomphLibError(
         err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
     }
  
     // Finish using the solver
     return;
   } // End of solve_helper

References abs(), beta, oomph::DoubleVector::built(), conj(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), hessenberg(), i, Global_Variables::Iterations, j, k, oomph::BlackBoxFDNewtonSolver::Max_iter, oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, s, BiharmonicTestFunctions1::solution(), sqrt(), oomph::Global_string_for_annotation::string(), oomph::TimingHelpers::timer(), and w.

Referenced by oomph::HelmholtzGMRESMG< MATRIX >::resolve(), and oomph::HelmholtzGMRESMG< MATRIX >::solve().

◆ update()

template<typename MATRIX >

void oomph::HelmholtzGMRESMG< MATRIX >::update	(	const unsigned &	k,
		const Vector< Vector< std::complex< double >>> &	hessenberg,
		const Vector< std::complex< double >> &	s,
		const Vector< Vector< DoubleVector >> &	v,
		Vector< DoubleVector > &	x
	)

inlineprotected

Helper function to update the result vector.

     {
       // Make a local copy of s
       Vector<std::complex<double>> y(s);
  
       //-----------------------------------------------------------------
       // The Hessenberg matrix should be an upper triangular matrix at
       // this point (although from its storage it would appear to be a
       // lower triangular matrix since the indexing has been reversed)
       // so finding the minimiser of J(y)=min||s-R_m*y|| where R_m is
       // the matrix R in the QR factorisation of the Hessenberg matrix.
       // Therefore, to obtain y we simply need to use a backwards
       // substitution. Note: The implementation here may appear to be
       // somewhat confusing as the indexing in the Hessenberg matrix is
       // reversed. This implementation of a backwards substitution does
       // not run along the columns of the triangular matrix but rather
       // up the rows.
       //-----------------------------------------------------------------
  
       // The outer loop is a loop over the columns of the Hessenberg matrix
       // since the indexing is reversed
       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] /= hessenberg[i][i];
  
         // The inner loop is a loop over the rows of the Hessenberg matrix
         for (int j = i - 1; j >= 0; j--)
         {
           // Update the j-th entry of y
           y[j] -= hessenberg[i][j] * y[i];
         }
       } // for (int i=int(k);i>=0;i--)
  
       // Calculate the number of entries in x (simply use the real part as
       // both the real and imaginary part should have the same length)
       unsigned n_row = x[0].nrow();
  
       // Build a temporary vector with entries initialised to 0.0
       Vector<DoubleVector> block_z(2);
  
       // Build a temporary vector with entries initialised to 0.0
       Vector<DoubleVector> block_temp(2);
  
       // Build the distributions
       for (unsigned dof_type = 0; dof_type < 2; dof_type++)
       {
         // Build the (dof_type)-th vector
         block_z[dof_type].build(x[0].distribution_pt(), 0.0);
  
         // Build the (dof_type)-th vector
         block_temp[dof_type].build(x[0].distribution_pt(), 0.0);
       }
  
       // Get access to the underlying values
       double* block_temp_r_pt = block_temp[0].values_pt();
  
       // Get access to the underlying values
       double* block_temp_c_pt = block_temp[1].values_pt();
  
       // Calculate x=Vy
       for (unsigned j = 0; j <= k; j++)
       {
         // Get access to j-th column of Z_m
         const double* vj_r_pt = v[j][0].values_pt();
  
         // Get access to j-th column of Z_m
         const double* vj_c_pt = v[j][1].values_pt();
  
         // Loop over the entries in x and update them
         for (unsigned i = 0; i < n_row; i++)
         {
           // Update the real part of the i-th entry in x
           block_temp_r_pt[i] +=
             (vj_r_pt[i] * y[j].real()) - (vj_c_pt[i] * y[j].imag());
  
           // Update the imaginary part of the i-th entry in x
           block_temp_c_pt[i] +=
             (vj_c_pt[i] * y[j].real()) + (vj_r_pt[i] * y[j].imag());
         }
       } // 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
         for (unsigned dof_type = 0; dof_type < 2; dof_type++)
         {
           // Update
           x[dof_type] += block_temp[dof_type];
         }
       }
       // If we're using RHS preconditioning
       else
       {
         // Create a temporary vector
         DoubleVector temp(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Copy block vectors block_temp back to temp
         this->return_block_vectors(block_temp, temp);
  
         // Create a temporary vector
         DoubleVector z(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Copy block vectors block_z back to z
         this->return_block_vectors(block_z, z);
  
         // Since we're using RHS preconditioning the preconditioner is applied
         // to the solution vector
         preconditioner_pt()->preconditioner_solve(temp, z);
  
         // Split up the solution vector into DoubleVectors, whose entries are
         // arranged to match the matrix blocks and assign it
         this->get_block_vectors(z, block_z);
  
         // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
         // sparse linear systems", p.284]
         for (unsigned dof_type = 0; dof_type < 2; dof_type++)
         {
           // Update
           x[dof_type] += block_z[dof_type];
         }
       } // if(Preconditioner_LHS) else
     } // End of update

protected

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

Referenced by oomph::HelmholtzGMRESMG< MATRIX >::resolve(), oomph::HelmholtzGMRESMG< MATRIX >::solve(), and oomph::HelmholtzFGMRESMG< MATRIX >::solve().

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

complex_smoother.h

Public Member Functions

Protected Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

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

Constructor & Destructor Documentation

◆ HelmholtzGMRESMG() [1/2]

◆ ~HelmholtzGMRESMG()

◆ HelmholtzGMRESMG() [2/2]

Member Function Documentation

◆ apply_plane_rotation()

◆ clean_up_memory()

◆ complex_matrix_multiplication()

◆ disable_resolve()

◆ generate_plane_rotation()

◆ iterations()

◆ operator=()

◆ preconditioner_solve()

◆ resolve()

◆ set_preconditioner_LHS()

◆ set_preconditioner_RHS()

◆ setup()

◆ solve() [1/3]

◆ solve() [2/3]

◆ solve() [3/3]

◆ solve_helper()

◆ update()

Member Data Documentation

◆ Iterations

◆ Matrices_storage_pt

◆ Matrix_can_be_deleted

◆ Preconditioner_LHS

◆ Resolving

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