The conjugate gradient method. More...

#include <iterative_linear_solver.h>

Inheritance diagram for oomph::BiCGStab< MATRIX >:

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

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

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

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

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

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

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

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

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

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

Public Member Functions inherited from oomph::IterativeLinearSolver
	IterativeLinearSolver ()

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

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

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

Preconditioner *&	preconditioner_pt ()
	Access function to preconditioner. More...

Preconditioner *const &	preconditioner_pt () const
	Access function to preconditioner (const version) More...

double &	tolerance ()
	Access to convergence tolerance. More...

unsigned &	max_iter ()
	Access to max. number of iterations. More...

void	enable_doc_convergence_history ()
	Enable documentation of the convergence history. More...

void	disable_doc_convergence_history ()
	Disable documentation of the convergence history. More...

void	open_convergence_history_file_stream (const std::string &file_name, const std::string &zone_title="")

void	close_convergence_history_file_stream ()
	Close convergence history output stream. More...

double	jacobian_setup_time () const

double	linear_solver_solution_time () const
	return the time taken to solve the linear system More...

virtual double	preconditioner_setup_time () const
	returns the the time taken to setup the preconditioner More...

void	enable_setup_preconditioner_before_solve ()
	Setup the preconditioner before the solve. More...

void	disable_setup_preconditioner_before_solve ()
	Don't set up the preconditioner before the solve. More...

void	enable_error_after_max_iter ()
	Throw an error if we don't converge within max_iter. More...

void	disable_error_after_max_iter ()
	Don't throw an error if we don't converge within max_iter (default). More...

void	enable_iterative_solver_as_preconditioner ()

void	disable_iterative_solver_as_preconditioner ()

Public Member Functions inherited from oomph::LinearSolver
	LinearSolver ()
	Empty constructor, initialise the member data. More...

	LinearSolver (const LinearSolver &dummy)=delete
	Broken copy constructor. More...

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

virtual	~LinearSolver ()
	Empty virtual destructor. More...

void	enable_doc_time ()
	Enable documentation of solve times. More...

void	disable_doc_time ()
	Disable documentation of solve times. More...

bool	is_doc_time_enabled () const
	Is documentation of solve times enabled? More...

bool	is_resolve_enabled () const
	Boolean flag indicating if resolves are enabled. More...

virtual void	enable_resolve ()

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

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

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

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

virtual void	enable_computation_of_gradient ()

void	disable_computation_of_gradient ()

void	reset_gradient ()

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

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

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

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

virtual	~DistributableLinearAlgebraObject ()
	Destructor. More...

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

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

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

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

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

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

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

bool	distribution_built () const

void	build_distribution (const LinearAlgebraDistribution *const dist_pt)

void	build_distribution (const LinearAlgebraDistribution &dist)

Private Member Functions
void	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...

Private Attributes
unsigned	Iterations
	Number of iterations taken. More...

MATRIX *	Matrix_pt
	Pointer to 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::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::BiCGStab< MATRIX >

The conjugate gradient method.

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

Constructor & Destructor Documentation

◆ BiCGStab() [1/2]

template<typename MATRIX >

oomph::BiCGStab< MATRIX >::BiCGStab ( )

inline

Constructor.

       : Iterations(0),
         Matrix_pt(0),
         Resolving(false),
         Matrix_can_be_deleted(true)
     {
     }

◆ ~BiCGStab()

template<typename MATRIX >

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

inlinevirtual

Destructor (cleanup storage)

     {
       clean_up_memory();
     }

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

◆ BiCGStab() [2/2]

template<typename MATRIX >

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

delete

Broken copy constructor.

Member Function Documentation

◆ clean_up_memory()

template<typename MATRIX >

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

inlineprivatevirtual

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

Reimplemented from oomph::LinearSolver.

     {
       if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
       {
         delete Matrix_pt;
         Matrix_pt = 0;
       }
     }

References oomph::BiCGStab< MATRIX >::Matrix_can_be_deleted, and oomph::BiCGStab< MATRIX >::Matrix_pt.

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

◆ disable_resolve()

template<typename MATRIX >

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

inlinevirtual

Overload disable resolve so that it cleans up memory too.

Reimplemented from oomph::LinearSolver.

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

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

◆ iterations()

template<typename MATRIX >

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

inlinevirtual

Number of iterations taken.

Implements oomph::IterativeLinearSolver.

     {
       return Iterations;
     }

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

◆ operator=()

template<typename MATRIX >

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

delete

Broken assignment operator.

◆ resolve()

template<typename MATRIX >

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

virtual

Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here. Solution is returned in the vector result.

//////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here. Solution is returned in the vector result.

Reimplemented from oomph::LinearSolver.

   {
     // We are re-solving
     Resolving = true;
  
 #ifdef PARANOID
     if (Matrix_pt == 0)
     {
       throw OomphLibError("No matrix was stored -- cannot re-solve",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Call linear algebra-style solver
     this->solve(Matrix_pt, rhs, result);
  
     // Reset re-solving flag
     Resolving = false;
   }

References OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and solve.

◆ solve() [1/3]

template<typename MATRIX >

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

inlinevirtual

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

Reimplemented from oomph::LinearSolver.

     {
       // Store the matrix if required
       if ((Enable_resolve) && (!Resolving))
       {
         Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
  
         // Matrix has been passed in from the outside so we must not
         // delete it
         Matrix_can_be_deleted = false;
       }
  
       // set the distribution
       if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt))
       {
         // the solver has the same distribution as the matrix if possible
         this->build_distribution(
           dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
             ->distribution_pt());
       }
       else
       {
         // the solver has the same distribution as the RHS
         this->build_distribution(rhs.distribution_pt());
       }
  
       // Call the helper function
       this->solve_helper(matrix_pt, rhs, solution);
     }

References oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Enable_resolve, oomph::BiCGStab< MATRIX >::Matrix_can_be_deleted, oomph::BiCGStab< MATRIX >::Matrix_pt, oomph::BiCGStab< MATRIX >::Resolving, BiharmonicTestFunctions1::solution(), and oomph::BiCGStab< MATRIX >::solve_helper().

◆ solve() [2/3]

template<typename MATRIX >

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

inlinevirtual

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

Reimplemented from oomph::LinearSolver.

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

References oomph::LinearSolver::solve().

◆ solve() [3/3]

template<typename MATRIX >

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

virtual

Solver: Takes pointer to problem and returns the results vector which contains the solution of the linear system defined by the problem's fully assembled Jacobian and residual vector.

Implements oomph::LinearSolver.

   {
     // Initialise timer
 #ifdef OOMPH_HAS_MPI
     double t_start = TimingHelpers::timer();
 #else
     clock_t t_start = clock();
 #endif
  
     // We're not re-solving
     Resolving = false;
  
     // Get rid of any previously stored data
     clean_up_memory();
  
     // Get Jacobian matrix in format specified by template parameter
     // and nonlinear residual vector
     Matrix_pt = new MATRIX;
     DoubleVector f;
     problem_pt->get_jacobian(f, *Matrix_pt);
  
  
     // We've made the matrix, we can delete it...
     Matrix_can_be_deleted = true;
  
     // Doc time for setup
 #ifdef OOMPH_HAS_MPI
     double t_end = TimingHelpers::timer();
     Jacobian_setup_time = t_end - t_start;
 #else
     clock_t t_end = clock();
     Jacobian_setup_time = double(t_end - t_start) / CLOCKS_PER_SEC;
 #endif
  
     if (Doc_time)
     {
       oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                  << std::endl;
     }
  
     // set the distribution
     if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt))
     {
       // the solver has the same distribution as the matrix if possible
       this->build_distribution(
         dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt)
           ->distribution_pt());
     }
     else
     {
       // the solver has the same distribution as the RHS
       this->build_distribution(f.distribution_pt());
     }
  
     // Call linear algebra-style solver
     if ((result.built()) &&
         (!(*result.distribution_pt() == *this->distribution_pt())))
     {
       LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
       result.build(this->distribution_pt(), 0.0);
       this->solve_helper(Matrix_pt, f, result);
       result.redistribute(&temp_global_dist);
     }
     else
     {
       this->solve_helper(Matrix_pt, f, result);
     }
  
     // Kill matrix unless it's still required for resolve
     if (!Enable_resolve) clean_up_memory();
   };

References oomph::DoubleVector::build(), oomph::DoubleVector::built(), oomph::TerminateHelper::clean_up_memory(), oomph::DistributableLinearAlgebraObject::distribution_pt(), f(), oomph::Problem::get_jacobian(), oomph::oomph_info, oomph::DoubleVector::redistribute(), and oomph::TimingHelpers::timer().

◆ solve_helper()

template<typename MATRIX >

void oomph::BiCGStab< 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. Algorithm and variable names based on "Numerical Linear Algebra for High-Performance Computers" by Dongarra, Duff, Sorensen & van der Vorst. SIAM (1998), page 185.

   {
 #ifdef PARANOID
     // check that the rhs vector is setup
     if (!rhs.built())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The vectors rhs must be setup";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
     // check that the matrix is square
     if (matrix_pt->nrow() != matrix_pt->ncol())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The matrix at matrix_pt must be square.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
     // check that the matrix and the rhs vector have the same nrow()
     if (matrix_pt->nrow() != rhs.nrow())
     {
       std::ostringstream error_message_stream;
       error_message_stream
         << "The matrix and the rhs vector must have the same number of rows.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
     // if the matrix is distributable then it too should have the same
     // communicator as the rhs vector
     DistributableLinearAlgebraObject* dist_matrix_pt =
       dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
     if (dist_matrix_pt != 0)
     {
       if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "The matrix matrix_pt must have the same communicator as the "
              "vectors"
           << " rhs and result must have the same communicator";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
     // if the matrix is not distributable then it the rhs vector should not be
     // distributed
     else
     {
       if (rhs.distribution_pt()->distributed())
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "The matrix (matrix_pt) is not distributable and therefore the rhs"
           << " vector must not be distributed";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
     // if the result vector is setup then check it has the same distribution
     // as the rhs
     if (solution.built())
     {
       if (!(*solution.distribution_pt() == *rhs.distribution_pt()))
       {
         std::ostringstream error_message_stream;
         error_message_stream << "The solution vector distribution has been "
                                 "setup; it must have the "
                              << "same distribution as the rhs vector.";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // setup the solution if it is not
     if (!solution.distribution_pt()->built())
     {
       solution.build(this->distribution_pt(), 0.0);
     }
     // zero
     else
     {
       solution.initialise(0.0);
     }
  
     // Get number of dofs
     // unsigned n_dof=rhs.size();
     unsigned nrow_local = this->nrow_local();
  
     // Time solver
 #ifdef OOMPH_HAS_MPI
     double t_start = TimingHelpers::timer();
 #else
     clock_t t_start = clock();
 #endif
  
     // Initialise: Zero initial guess so the initial residual is
     // equal to the RHS, i.e. the nonlinear residual
     DoubleVector residual(rhs);
     double residual_norm = residual.norm();
     double rhs_norm = residual_norm;
     if (rhs_norm == 0.0) rhs_norm = 1.0;
     DoubleVector x(rhs.distribution_pt(), 0.0);
  
     // Hat residual by copy operation
     DoubleVector r_hat(residual);
  
     // Normalised residual
     double normalised_residual_norm = residual_norm / rhs_norm;
  
     // if required will document convergence history to screen or file (if
     // stream open)
     if (Doc_convergence_history)
     {
       if (!Output_file_stream.is_open())
       {
         oomph_info << 0 << " " << normalised_residual_norm << std::endl;
       }
       else
       {
         Output_file_stream << 0 << " " << normalised_residual_norm << std::endl;
       }
     }
  
     // Check immediate convergence
     if (normalised_residual_norm < Tolerance)
     {
       if (Doc_time)
       {
         oomph_info << "BiCGStab converged immediately" << std::endl;
       }
       solution = x;
  
       // Doc time for solver
       double t_end = TimingHelpers::timer();
       Solution_time = t_end - t_start;
  
       if (Doc_time)
       {
         oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                    << std::endl;
       }
       return;
     }
  
     // Setup preconditioner only if we're not re-solving
     if (!Resolving)
     {
       // only setup the preconditioner if required
       if (Setup_preconditioner_before_solve)
       {
         // Setup preconditioner from the Jacobian matrix
         double t_start_prec = TimingHelpers::timer();
  
         preconditioner_pt()->setup(matrix_pt);
  
         // Doc time for setup of preconditioner
         double t_end_prec = TimingHelpers::timer();
         Preconditioner_setup_time = t_end_prec - t_start_prec;
  
         if (Doc_time)
         {
           oomph_info << "Time for setup of preconditioner  [sec]: "
                      << Preconditioner_setup_time << std::endl;
         }
       }
     }
     else
     {
       if (Doc_time)
       {
         oomph_info << "Setup of preconditioner is bypassed in resolve mode"
                    << std::endl;
       }
     }
  
     // Some auxiliary variables
     double rho_prev = 0.0;
     double alpha = 0.0;
     double omega = 0.0;
     double rho, beta, dot_prod, dot_prod_tt, dot_prod_ts;
     double s_norm, r_norm;
  
     // Some vectors
     DoubleVector p(this->distribution_pt(), 0.0),
       p_hat(this->distribution_pt(), 0.0), v(this->distribution_pt(), 0.0),
       z(this->distribution_pt(), 0.0), t(this->distribution_pt(), 0.0),
       s(this->distribution_pt(), 0.0);
  
     // Loop over max. number of iterations
     for (unsigned iter = 1; iter <= Max_iter; iter++)
     {
       // Dot product for rho
       rho = r_hat.dot(residual);
  
       // Breakdown?
       if (rho == 0.0)
       {
         oomph_info << "BiCGStab has broken down after " << iter << " iterations"
                    << std::endl;
         oomph_info << "Returning with current normalised residual of "
                    << normalised_residual_norm << std::endl;
       }
  
       // First step is different
       if (iter == 1)
       {
         p = residual;
       }
       else
       {
         beta = (rho / rho_prev) * (alpha / omega);
         for (unsigned i = 0; i < nrow_local; i++)
         {
           p[i] = residual[i] + beta * (p[i] - omega * v[i]);
         }
       }
  
       // Apply precondtitioner: p_hat=P^-1*p
       preconditioner_pt()->preconditioner_solve(p, p_hat);
  
       // Matrix vector product: v=A*p_hat
       matrix_pt->multiply(p_hat, v);
       dot_prod = r_hat.dot(v);
       alpha = rho / dot_prod;
       for (unsigned i = 0; i < nrow_local; i++)
       {
         s[i] = residual[i] - alpha * v[i];
       }
       s_norm = s.norm();
  
       // Normalised residual
       normalised_residual_norm = s_norm / rhs_norm;
  
       // if required will document convergence history to screen or file (if
       // stream open)
       if (Doc_convergence_history)
       {
         if (!Output_file_stream.is_open())
         {
           oomph_info << double(iter - 0.5) << " " << normalised_residual_norm
                      << std::endl;
         }
         else
         {
           Output_file_stream << double(iter - 0.5) << " "
                              << normalised_residual_norm << std::endl;
         }
       }
  
       // Converged?
       if (normalised_residual_norm < Tolerance)
       {
         for (unsigned i = 0; i < nrow_local; i++)
         {
           solution[i] = x[i] + alpha * p_hat[i];
         }
  
         if (Doc_time)
         {
           oomph_info << std::endl;
           oomph_info << "BiCGStab converged. Normalised residual norm: "
                      << normalised_residual_norm << std::endl;
           oomph_info << "Number of iterations to convergence: " << iter
                      << std::endl;
           oomph_info << std::endl;
         }
  
         // Store number of iterations taken
         Iterations = iter;
  
         // Doc time for solver
         double t_end = TimingHelpers::timer();
         Solution_time = t_end - t_start;
  
         if (Doc_time)
         {
           oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                      << std::endl;
         }
  
         return;
       }
  
       // Apply precondtitioner: z=P^-1*s
       preconditioner_pt()->preconditioner_solve(s, z);
  
       // Matrix vector product: t=A*z
       matrix_pt->multiply(z, t);
       dot_prod_ts = t.dot(s);
       dot_prod_tt = t.dot(t);
       omega = dot_prod_ts / dot_prod_tt;
       for (unsigned i = 0; i < nrow_local; i++)
       {
         x[i] += alpha * p_hat[i] + omega * z[i];
         residual[i] = s[i] - omega * t[i];
       }
       r_norm = residual.norm();
       rho_prev = rho;
  
       // Check convergence again
       normalised_residual_norm = r_norm / rhs_norm;
  
       // if required will document convergence history to screen or file (if
       // stream open)
       if (Doc_convergence_history)
       {
         if (!Output_file_stream.is_open())
         {
           oomph_info << iter << " " << normalised_residual_norm << std::endl;
         }
         else
         {
           Output_file_stream << iter << " " << normalised_residual_norm
                              << std::endl;
         }
       }
  
  
       if (normalised_residual_norm < Tolerance)
       {
         if (Doc_time)
         {
           oomph_info << std::endl;
           oomph_info << "BiCGStab converged. Normalised residual norm: "
                      << normalised_residual_norm << std::endl;
           oomph_info << "Number of iterations to convergence: " << iter
                      << std::endl;
           oomph_info << std::endl;
         }
         solution = x;
  
         // Store the number of itertions taken.
         Iterations = iter;
  
         // Doc time for solver
         double t_end = TimingHelpers::timer();
         Solution_time = t_end - t_start;
  
         if (Doc_time)
         {
           oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                      << std::endl;
         }
         return;
       }
  
  
       // Breakdown: Omega has to be >0 for to be able to continue
       if (omega == 0.0)
       {
         oomph_info << std::endl;
         oomph_info << "BiCGStab breakdown with omega=0.0. "
                    << "Normalised residual norm: " << normalised_residual_norm
                    << std::endl;
         oomph_info << "Number of iterations so far: " << iter << std::endl;
         oomph_info << std::endl;
         solution = x;
  
         // Store the number of itertions taken.
         Iterations = iter;
  
         // Doc time for solver
         double t_end = TimingHelpers::timer();
         Solution_time = t_end - t_start;
  
         if (Doc_time)
         {
           oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                      << std::endl;
         }
         return;
       }
  
  
     } // end of iteration loop
  
  
     // No convergence
     oomph_info << std::endl;
     oomph_info << "BiCGStab did not converge to required tolerance! "
                << std::endl;
     oomph_info << "Returning with normalised residual norm: "
                << normalised_residual_norm << std::endl;
     oomph_info << "after " << Max_iter << " iterations." << std::endl;
     oomph_info << std::endl;
  
     solution = x;
  
     // Doc time for solver
     double t_end = TimingHelpers::timer();
     Solution_time = t_end - t_start;
  
     if (Doc_time)
     {
       oomph_info << "Time for solve with BiCGStab  [sec]: " << Solution_time
                  << std::endl;
     }
  
     if (Throw_error_after_max_iter)
     {
       std::string err = "Solver failed to converge and you requested an error";
       err += " on convergence failures.";
       throw OomphLibError(
         err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
     }
   }

References alpha, beta, oomph::DoubleVector::built(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::DoubleVector::dot(), i, Global_Variables::Iterations, oomph::BlackBoxFDNewtonSolver::Max_iter, oomph::DoubleMatrixBase::multiply(), oomph::DoubleMatrixBase::ncol(), oomph::DoubleVector::norm(), oomph::DistributableLinearAlgebraObject::nrow(), oomph::DoubleMatrixBase::nrow(), Eigen::internal::omega(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, p, s, BiharmonicTestFunctions1::solution(), oomph::Global_string_for_annotation::string(), plotPSD::t, oomph::TimingHelpers::timer(), v, and plotDoE::x.

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

Member Data Documentation

◆ Iterations

template<typename MATRIX >

unsigned oomph::BiCGStab< MATRIX >::Iterations

private

Number of iterations taken.

Referenced by oomph::BiCGStab< MATRIX >::iterations().

◆ Matrix_can_be_deleted

template<typename MATRIX >

bool oomph::BiCGStab< MATRIX >::Matrix_can_be_deleted

private

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

Referenced by oomph::BiCGStab< MATRIX >::clean_up_memory(), and oomph::BiCGStab< MATRIX >::solve().

◆ Matrix_pt

template<typename MATRIX >

MATRIX* oomph::BiCGStab< MATRIX >::Matrix_pt

private

Pointer to matrix.

Referenced by oomph::BiCGStab< MATRIX >::clean_up_memory(), and oomph::BiCGStab< MATRIX >::solve().

◆ Resolving

template<typename MATRIX >

bool oomph::BiCGStab< MATRIX >::Resolving

private

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

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

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

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

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

Constructor & Destructor Documentation

◆ BiCGStab() [1/2]

◆ ~BiCGStab()

◆ BiCGStab() [2/2]

Member Function Documentation

◆ clean_up_memory()

◆ disable_resolve()

◆ iterations()

◆ operator=()

◆ resolve()

◆ solve() [1/3]

◆ solve() [2/3]

◆ solve() [3/3]

◆ solve_helper()

Member Data Documentation

◆ Iterations

◆ Matrix_can_be_deleted

◆ Matrix_pt

◆ Resolving

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