#include <complex_smoother.h>

Inheritance diagram for oomph::HelmholtzFGMRESMG< MATRIX >:

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

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

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

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

void	set_preconditioner_LHS ()

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

Public Member Functions inherited from oomph::HelmholtzGMRESMG< MATRIX >
	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 (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)

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

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

Additional Inherited Members
Protected Member Functions inherited from oomph::HelmholtzGMRESMG< MATRIX >
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 inherited from oomph::HelmholtzGMRESMG< MATRIX >
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...

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

Detailed Description

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

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// The FGMRES method, i.e. the flexible variant of the GMRES method which allows for nonconstant preconditioners [see Saad Y, "Iterative methods for sparse linear systems", p.287]. Note, FGMRES can only cater to right preconditioning; if the user tries to switch to left preconditioning they will be notified of this

Constructor & Destructor Documentation

◆ HelmholtzFGMRESMG() [1/2]

template<typename MATRIX >

oomph::HelmholtzFGMRESMG< MATRIX >::HelmholtzFGMRESMG ( )

inline

Constructor (empty)

                         : HelmholtzGMRESMG<MATRIX>()
     {
       // Can only use RHS preconditioning
       this->Preconditioner_LHS = false;
     };

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

◆ ~HelmholtzFGMRESMG()

template<typename MATRIX >

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

inlinevirtual

Destructor (cleanup storage)

     {
       // Call the clean up function in the base class GMRES
       this->clean_up_memory();
     }

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

◆ HelmholtzFGMRESMG() [2/2]

template<typename MATRIX >

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

delete

Broken copy constructor.

Member Function Documentation

◆ operator=()

template<typename MATRIX >

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

delete

Broken assignment operator.

◆ set_preconditioner_LHS()

template<typename MATRIX >

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

inline

Overloaded function to let the user know that left preconditioning is not possible with FGMRES, only right preconditioning

     {
       // Create an output stream
       std::ostringstream error_message_stream;
  
       // Create an error message
       error_message_stream << "FGMRES cannot use left preconditioning. It is "
                            << "only capable of using right preconditioning."
                            << std::endl;
  
       // Throw the error message
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     } // End of set_preconditioner_LHS

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ solve()

template<typename MATRIX >

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

Reimplemented from oomph::HelmholtzGMRESMG< 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
       this->Resolving = false;
  
       // Get rid of any previously stored data
       this->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...
       this->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
       this->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
         this->Matrices_storage_pt[i] = new CRDoubleMatrix;
  
         // Extract the required blocks, i.e. the first column
         this->get_block(i, j, *(this->Matrices_storage_pt[i]));
       }
  
       // Doc time for setup
       double t_end = TimingHelpers::timer();
       this->Jacobian_setup_time = t_end - t_start;
  
       if (this->Doc_time)
       {
         oomph_info << "\nTime for setup of block Jacobian [sec]: "
                    << this->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 (!(this->Enable_resolve))
       {
         // Clean up anything left in memory
         this->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::HelmholtzFGMRESMG< MATRIX >::solve_helper(), oomph::TimingHelpers::timer(), and tmp.

◆ solve_helper()

template<typename MATRIX >

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

private

General interface to solve function.

/////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the linear system. based on the algorithm presented in Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, Barrett, Berry et al, SIAM, 2006 and the implementation in the IML++ library : http://math.nist.gov/iml++/

   {
     // 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 = this->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 (this->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: "
                            << this->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*>(
             this->Matrices_storage_pt[dof_type]) != 0)
       {
         if (dynamic_cast<DistributableLinearAlgebraObject*>(
               this->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 (!(this->Resolving))
     {
       // Only set up the preconditioner before solve if required
       if (this->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
         this->preconditioner_pt()->setup(dynamic_cast<MATRIX*>(matrix_pt));
  
         // Doc time for setup of preconditioner
         double t_end_prec = TimingHelpers::timer();
         this->Preconditioner_setup_time = t_end_prec - t_start_prec;
  
         // If time documentation is enabled
         if (this->Doc_time)
         {
           // Output the time taken
           oomph_info << "Time for setup of preconditioner [sec]: "
                      << this->Preconditioner_setup_time << std::endl;
         }
       }
     }
     else
     {
       // If time documentation is enabled
       if (this->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);
     }
  
     // 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];
     }
  
     // 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 (this->Doc_convergence_history)
     {
       // If an output file which is open isn't provided then output to screen
       if (!(this->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
         this->Output_file_stream << 0 << " " << resid << std::endl;
       }
     } // if (Doc_convergence_history)
  
     // If the GMRES algorithm converges immediately
     if (resid <= this->Tolerance)
     {
       // If time documentation is enabled
       if (this->Doc_time)
       {
         // Notify the user
         oomph_info << "FGMRES 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);
  
     // Create a vector of DoubleVectors (stores the preconditioned vectors)
     Vector<Vector<DoubleVector>> block_z;
  
     // Resize the number of vectors needed
     block_z.resize(n_row + 1);
  
     // Resize each Vector of DoubleVectors to store the real and imaginary
     // part of a given vector
     for (unsigned i = 0; i < n_row + 1; i++)
     {
       // Create space for two DoubleVector objects in each Vector
       block_v[i].resize(n_dof_types);
  
       // Create space for two DoubleVector objects in each Vector
       block_z[i].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* block_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
         block_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 < this->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_j=Jz_j where z_j=M^{-1}v_j (RHS prec.)
       {
         // Create a temporary vector
         DoubleVector vj(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Create a temporary vector
         DoubleVector zj(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Create a temporary vector
         DoubleVector w(IterativeLinearSolver::distribution_pt(), 0.0);
  
         // Create two DoubleVectors
         for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
         {
           // Build the k-th part of the j-th preconditioning result vector
           block_z[j][dof_type].build(this->block_distribution_pt(dof_type),
                                      0.0);
         }
  
         // Copy the real and imaginary part of v[j] into one vector, vj
         this->return_block_vectors(block_v[j], vj);
  
         // Calculate z_j=M^{-1}v_j by saad p270
         this->preconditioner_pt()->preconditioner_solve(vj, zj);
  
         // Copy zj into z[j][0] and z[j][1]
         this->get_block_vectors(zj, block_z[j]);
  
         // Calculate w_j=J*(M^{-1}v_j). Note, we cannot use inbuilt complex
         // matrix algebra here as we're using distributed vectors
         this->complex_matrix_multiplication(
           this->Matrices_storage_pt, block_z[j], block_w);
       } // Calculate 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* block_vi_r_pt = block_v[i][0].values_pt();
  
         // For fast access
         const double* block_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(block_vi_r_pt[k], block_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() * block_vi_r_pt[k] -
                               hessenberg[j][i].imag() * block_vi_c_pt[k]);
  
           // Update the imaginary part of the k-th entry of w
           block_w_c_pt[k] -= (hessenberg[j][i].real() * block_vi_c_pt[k] +
                               hessenberg[j][i].imag() * block_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* block_v_r_pt = block_v[j + 1][0].values_pt();
  
         // For fast access
         double* block_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
           block_v_r_pt[k] = block_w_r_pt[k] / h_subdiag_val;
  
           // Similarly, the i-th entry of the imaginary component is given by
           block_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)
         this->apply_plane_rotation(
           hessenberg[j][k], hessenberg[j][k + 1], cs[k], sn[k]);
       }
  
       // Now calculate the entries of the latest Givens rotation matrix
       this->generate_plane_rotation(
         hessenberg[j][j], hessenberg[j][j + 1], cs[j], sn[j]);
  
       // Apply the plane rotation using the newly calculated entries
       this->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||
       this->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 (this->Doc_convergence_history)
       {
         // If an output file which is open isn't provided then output to screen
         if (!(this->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
           this->Output_file_stream << j + 1 << " " << resid << std::endl;
         }
       } // if (Doc_convergence_history)
  
       // If the required tolerance has been met
       if (resid < this->Tolerance)
       {
         // Store the number of iterations taken
         this->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_z, block_solution);
  
         // Copy the vectors in block_solution to solution
         this->return_block_vectors(block_solution, solution);
  
         // If time documentation was enabled
         if (this->Doc_time)
         {
           // Output the current normalised residual norm
           oomph_info << "\nFGMRES 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
         this->Solution_time = t_end - t_start;
  
         // If time documentation was enabled
         if (this->Doc_time)
         {
           // Output the time taken to solve the problem using GMRES
           oomph_info << "Time for solve with FGMRES [sec]: "
                      << this->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
     this->Iterations = this->Max_iter;
  
     // Only update if we actually did something
     if (this->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(this->Max_iter - 1, hessenberg, s, block_z, block_solution);
  
       // Copy the vectors in block_solution to solution
       this->return_block_vectors(block_solution, solution);
     }
  
     // 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 < this->Tolerance)
     {
       // If time documentation is enabled
       if (this->Doc_time)
       {
         // Notify the user
         oomph_info << "\nFGMRES converged (2). Normalised residual norm: "
                    << resid << "\nNumber of iterations to convergence: "
                    << this->Iterations << "\n"
                    << std::endl;
       }
  
       // End the timer
       double t_end = TimingHelpers::timer();
  
       // Calculate the time taken for the solver
       this->Solution_time = t_end - t_start;
  
       // If time documentation is enabled
       if (this->Doc_time)
       {
         oomph_info << "Time for solve with FGMRES [sec]: "
                    << this->Solution_time << std::endl;
       }
       return;
     }
  
     // Otherwise GMRES failed convergence
     oomph_info << "\nFGMRES did not converge to required tolerance! "
                << "\nReturning with normalised residual norm: " << resid
                << "\nafter " << this->Max_iter << " iterations.\n"
                << std::endl;
  
     // Throw an error if requested
     if (this->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(), s, BiharmonicTestFunctions1::solution(), sqrt(), oomph::Global_string_for_annotation::string(), oomph::TimingHelpers::timer(), and w.

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

◆ update()

template<typename MATRIX >

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

inlineprivate

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_update(2);
  
       // Build the distributions
       for (unsigned dof_type = 0; dof_type < 2; dof_type++)
       {
         // Build the (dof_type)-th vector of block_update
         block_update[dof_type].build(x[0].distribution_pt(), 0.0);
       }
  
       // Get access to the underlying values
       double* block_update_r_pt = block_update[0].values_pt();
  
       // Get access to the underlying values
       double* block_update_c_pt = block_update[1].values_pt();
  
       // Calculate x=Vy
       for (unsigned j = 0; j <= k; j++)
       {
         // Get access to j-th column of Z_m
         const double* z_mj_r_pt = z_m[j][0].values_pt();
  
         // Get access to j-th column of Z_m
         const double* z_mj_c_pt = z_m[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_update_r_pt[i] +=
             (z_mj_r_pt[i] * y[j].real()) - (z_mj_c_pt[i] * y[j].imag());
  
           // Update the imaginary part of the i-th entry in x
           block_update_c_pt[i] +=
             (z_mj_c_pt[i] * y[j].real()) + (z_mj_r_pt[i] * y[j].imag());
         }
       } // for (unsigned j=0;j<=k;j++)
  
       // 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 the solution
         x[dof_type] += block_update[dof_type];
       }
     } // End of update