The GMRES method rewritten for complex matrices. More...

#include <complex_smoother.h>

Inheritance diagram for oomph::ComplexGMRES< MATRIX >:

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

	~ComplexGMRES ()
	Empty destructor. More...

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

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

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

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

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

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

void	complex_smoother_setup (Vector< CRDoubleMatrix * > helmholtz_matrix_pt)
	Setup: Pass pointer to the matrix and store in cast form. More...

void	complex_smoother_solve (const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)

Public Member Functions inherited from oomph::HelmholtzSmoother
	HelmholtzSmoother ()
	Empty constructor. More...

virtual	~HelmholtzSmoother ()
	Virtual empty destructor. More...

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

template<typename MATRIX >
void	check_validity_of_solve_helper_inputs (CRDoubleMatrix const &real_matrix_pt, CRDoubleMatrix const &imag_matrix_pt, const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution, const double &n_dof)

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 (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)

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

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

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

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

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

virtual void	enable_computation_of_gradient ()

void	disable_computation_of_gradient ()

void	reset_gradient ()

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

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

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

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

virtual	~DistributableLinearAlgebraObject ()
	Destructor. More...

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

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

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

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

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

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

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

bool	distribution_built () const

void	build_distribution (const LinearAlgebraDistribution *const dist_pt)

void	build_distribution (const LinearAlgebraDistribution &dist)

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

void	complex_solve_helper (const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
	This is where the actual work is done. More...

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)

Private 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

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

Protected Attributes inherited from oomph::HelmholtzSmoother
bool	Use_as_smoother

Protected Attributes inherited from oomph::IterativeLinearSolver
bool	Doc_convergence_history

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

double	Tolerance
	Convergence tolerance. More...

unsigned	Max_iter
	Maximum number of iterations. More...

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

double	Jacobian_setup_time
	Jacobian setup time. More...

double	Solution_time
	linear solver solution time More...

double	Preconditioner_setup_time
	Preconditioner setup time. More...

bool	Setup_preconditioner_before_solve

bool	Throw_error_after_max_iter

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

bool	First_time_solve_when_used_as_preconditioner

Protected Attributes inherited from oomph::LinearSolver
bool	Enable_resolve

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

bool	Compute_gradient

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

DoubleVector	Gradient_for_glob_conv_newton_solve

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

Detailed Description

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

The GMRES method rewritten for complex matrices.

Constructor & Destructor Documentation

◆ ComplexGMRES() [1/2]

template<typename MATRIX >

oomph::ComplexGMRES< MATRIX >::ComplexGMRES ( )

inline

Constructor.

       : Iterations(0),
         Matrices_storage_pt(0),
         Resolving(false),
         Matrix_can_be_deleted(true)
     {
     } // End of ComplexGMRES constructor

◆ ~ComplexGMRES()

template<typename MATRIX >

oomph::ComplexGMRES< MATRIX >::~ComplexGMRES ( )

inline

Empty destructor.

     {
       // Run clean_up_memory
       clean_up_memory();
     } // End of ~ComplexGMRES

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

◆ ComplexGMRES() [2/2]

template<typename MATRIX >

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

delete

Broken copy constructor.

Member Function Documentation

◆ apply_plane_rotation()

template<typename MATRIX >

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

inlineprivate

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::ComplexGMRES< MATRIX >::clean_up_memory ( )

inlineprivatevirtual

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

Reimplemented from oomph::LinearSolver.

     {
       // 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::ComplexGMRES< MATRIX >::Matrices_storage_pt, and oomph::ComplexGMRES< MATRIX >::Matrix_can_be_deleted.

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

◆ complex_smoother_setup()

template<typename MATRIX >

void oomph::ComplexGMRES< MATRIX >::complex_smoother_setup ( Vector< CRDoubleMatrix * > helmholtz_matrix_pt )

inlinevirtual

Setup: Pass pointer to the matrix and store in cast form.

Implements oomph::HelmholtzSmoother.

     {
       // Assume the matrices have been passed in from the outside so we must
       // not delete it. This is needed to avoid pre- and post-smoothers
       // deleting the same matrix in the MG solver. If this was originally
       // set to TRUE then this will be sorted out in the other functions
       // from which this was called
       Matrix_can_be_deleted = false;
  
 #ifdef PARANOID
       // PARANOID check - Make sure the input has the right number of matrices
       if (helmholtz_matrix_pt.size() != 2)
       {
         std::ostringstream error_message_stream;
         error_message_stream << "Can only deal with two matrices. You have "
                              << helmholtz_matrix_pt.size() << " matrices."
                              << std::endl;
  
         // Throw an error
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       // Resize the storage for the system matrices
       Matrices_storage_pt.resize(2, 0);
  
       // Assign the real matrix pointer
       Matrices_storage_pt[0] = helmholtz_matrix_pt[0];
  
       // Assign the imaginary matrix pointers
       Matrices_storage_pt[1] = helmholtz_matrix_pt[1];
  
 #ifdef PARANOID
       // PARANOID check - Make sure that the constituent matrices have the same
       // number of rows
       if (Matrices_storage_pt[0]->nrow() != Matrices_storage_pt[1]->nrow())
       {
         std::ostringstream error_message_stream;
         error_message_stream << "Incompatible real and imag. matrix sizes.";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
       // PARANOID check - Make sure that the constituent matrices have the same
       // number of columns
       if (Matrices_storage_pt[0]->ncol() != Matrices_storage_pt[1]->ncol())
       {
         std::ostringstream error_message_stream;
         error_message_stream << "Incompatible real and imag. matrix sizes.";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
 #endif
     } // End of complex_smoother_setup

References oomph::ComplexGMRES< MATRIX >::Matrices_storage_pt, oomph::ComplexGMRES< MATRIX >::Matrix_can_be_deleted, oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ complex_smoother_solve()

template<typename MATRIX >

void oomph::ComplexGMRES< MATRIX >::complex_smoother_solve	(	const Vector< DoubleVector > &	rhs,
		Vector< DoubleVector > &	solution
	)

inlinevirtual

The smoother_solve function performs fixed number of iterations on the system A*result=rhs. The number of (smoothing) iterations is the same as the max. number of iterations in the underlying IterativeLinearSolver class.

Implements oomph::HelmholtzSmoother.

     {
       // If you use a smoother but you don't want to calculate the residual
       Use_as_smoother = true;
  
       // Use the helper function where the work is actually done
       complex_solve_helper(rhs, solution);
     } // End of complex_smoother_solve

References oomph::ComplexGMRES< MATRIX >::complex_solve_helper(), BiharmonicTestFunctions1::solution(), and oomph::HelmholtzSmoother::Use_as_smoother.

◆ complex_solve_helper()

template<typename MATRIX >

void oomph::ComplexGMRES< MATRIX >::complex_solve_helper	(	const Vector< DoubleVector > &	rhs,
		Vector< DoubleVector > &	solution
	)

private

This is where the actual work is done.

   {
     // Set the number of dof types (real and imaginary for this solver)
     unsigned n_dof_types = 2;
  
     // 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);
     }
  
     // Loop over the real and imaginary parts
     for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
     {
 #ifdef PARANOID
       // PARANOID check that if the matrix is distributable then it should not
       // be then it should not be distributed
       if (dynamic_cast<DistributableLinearAlgebraObject*>(
             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[dof_type].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[dof_type].nrow() != 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[dof_type].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[dof_type].built())
       {
         if (!(*rhs[dof_type].distribution_pt() ==
               *solution[dof_type].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[dof_type].built())
       {
         // Build the distribution
         solution[dof_type].build(this->distribution_pt(), 0.0);
       }
       // Otherwise initialise all entries to zero
       else
       {
         // Initialise the entries in the k-th vector in solution to zero
         solution[dof_type].initialise(0.0);
       }
     } // for (unsigned dof_type=0;dof_type<n_dof_types;dof_type++)
  
     // 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> 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
       w[dof_type].build(this->distribution_pt(), 0.0);
     }
  
     // 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> 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
       r[dof_type].build(this->distribution_pt(), 0.0);
     }
  
     // Store the value of b (the RHS vector) in r
     for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
     {
       // Build the distribution of the k-th constituent vector
       r[dof_type] = rhs[dof_type];
     }
  
     // Calculate the norm of the real part of r
     double norm_r = r[0].norm();
  
     // Calculate the norm of the imaginary part of r
     double norm_c = 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>> v;
  
     // Resize the number of vectors needed
     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
       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
       v[0][dof_type].build(this->distribution_pt(), 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 = v[0][dof_type].values_pt();
  
       // For fast access
       const double* r_pt = 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] = 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=J*v_j. Note, we cannot use inbuilt complex matrix algebra
       // here as we're using distributed vectors
       complex_matrix_multiplication(Matrices_storage_pt, v[j], w);
  
       // For fast access
       double* w_r_pt = w[0].values_pt();
  
       // For fast access
       double* w_c_pt = 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 = v[i][0].values_pt();
  
         // For fast access
         const double* vi_c_pt = 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(w_r_pt[k], w_c_pt[k]);
  
           // Update the value of H(i,j) noting we're computing a complex
           // inner product 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
           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
           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 = w[0].norm();
  
       // Calculate the norm of the imaginary part of w
       norm_c = 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 real part of the next orthogonal vector
       v[j + 1][0].build(this->distribution_pt(), 0.0);
  
       // Build the imaginary part of the next orthogonal vector
       v[j + 1][1].build(this->distribution_pt(), 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 = v[j + 1][0].values_pt();
  
         // For fast access
         double* v_c_pt = v[j + 1][1].values_pt();
  
         // For fast access
         const double* w_r_pt = w[0].values_pt();
  
         // For fast access
         const double* w_c_pt = 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] = w_r_pt[k] / h_subdiag_val;
  
           // Similarly, the i-th entry of the imaginary component is given by
           v_c_pt[k] = 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, v, 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, v, solution);
     }
  
     // 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;
     }
  
     // Finish using the solver
     return;
   } // End of complex_solve_helper

References abs(), beta, conj(), hessenberg(), i, Global_Variables::Iterations, j, k, oomph::BlackBoxFDNewtonSolver::Max_iter, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, s, BiharmonicTestFunctions1::solution(), sqrt(), oomph::TimingHelpers::timer(), v, and w.

Referenced by oomph::ComplexGMRES< MATRIX >::complex_smoother_solve().

◆ disable_resolve()

template<typename MATRIX >

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

inlinevirtual

Overload disable resolve so that it cleans up memory too.

Reimplemented from oomph::LinearSolver.

     {
       // Disable resolve using the LinearSolver function
       LinearSolver::disable_resolve();
  
       // Clean up anything kept in memory
       clean_up_memory();
     } // End of disable resolve

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

◆ generate_plane_rotation()

template<typename MATRIX >

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

inlineprivate

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} \cr -\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::ComplexGMRES< MATRIX >::iterations ( ) const

inlinevirtual

Number of iterations taken.

Implements oomph::IterativeLinearSolver.

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

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

◆ operator=()

template<typename MATRIX >

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

delete

Broken assignment operator.

◆ solve() [1/2]

template<typename MATRIX >

void oomph::ComplexGMRES< 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);
     } // End of solve

References oomph::LinearSolver::solve().

◆ solve() [2/2]

template<typename MATRIX >

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

     {
       // Write the error message into a string
       std::string error_message = "Solve function for class\n\n";
       error_message += "ComplexGMRES\n\n";
       error_message += "is deliberately broken to avoid the accidental \n";
       error_message += "use of the inappropriate C++ default. If you \n";
       error_message += "really need one for this class, write it yourself...\n";
  
       // Throw an error
       throw OomphLibError(
         error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     } // End of solve

References OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::Global_string_for_annotation::string().

◆ update()

template<typename MATRIX >

void oomph::ComplexGMRES< 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
	)

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();
  
       // We assume here that the vector x (which is passed in) is actually x_0
       // so we simply need to update its entries to calculate the solution, x
       // which is given by x=x_0+Vy.
       for (unsigned j = 0; j <= k; j++)
       {
         // For fast access (real part)
         const double* vj_r_pt = v[j][0].values_pt();
  
         // For fast access (imaginary part)
         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
           x[0][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
           x[1][i] += (vj_c_pt[i] * y[j].real()) + (vj_r_pt[i] * y[j].imag());
         }
       } // for (unsigned j=0;j<=k;j++)
     } // End of update