oomph::GS< CRDoubleMatrix > Class Reference

#include <iterative_linear_solver.h>

Inheritance diagram for oomph::GS< CRDoubleMatrix >:

Public Member Functions
	GS ()
	Constructor. More...
virtual	~GS ()
	Destructor (cleanup storage) More...
	GS (const GS &)=delete
	Broken copy constructor. More...
void	operator= (const GS &)=delete
	Broken assignment operator. More...
void	smoother_solve (const DoubleVector &rhs, DoubleVector &result)
void	disable_resolve ()
	Overload disable resolve so that it cleans up memory too. More...
void	smoother_setup (DoubleMatrixBase *matrix_pt)
	Set up the smoother for the matrix specified by the pointer. More...
void	setup_helper (DoubleMatrixBase *matrix_pt)
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)
double	preconditioner_setup_time () const
	Returns the time taken to set up the preconditioner. More...
unsigned	iterations () const
	Number of iterations taken. More...
Public Member Functions inherited from oomph::Smoother
	Smoother ()
	Empty constructor. More...
virtual	~Smoother ()
	Virtual empty destructor. More...
template<typename MATRIX >
void	check_validity_of_solve_helper_inputs (MATRIX *const &matrix_pt, const DoubleVector &rhs, 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...
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 ()
	Clean up data that's stored for resolve (if any has been stored) More...

Private Attributes
CRDoubleMatrix *	Matrix_pt
	System matrix pointer in the format specified by the template argument. More...
unsigned	Iterations
	Number of iterations taken. More...
bool	Resolving
bool	Matrix_can_be_deleted
Vector< int >	Index_of_diagonal_entries

Additional Inherited Members
Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void	clear_distribution ()
Protected Attributes inherited from oomph::Smoother
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

//////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////// Explicit template specialisation of the Gauss Seidel method for compressed row format matrices

Constructor & Destructor Documentation

GS() [1/2]

oomph::GS< CRDoubleMatrix >::GS ( )

inline

Constructor.

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

~GS()

virtual oomph::GS< CRDoubleMatrix >::~GS ( )

inlinevirtual

Destructor (cleanup storage)

    {
      clean_up_memory();
    }

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

GS() [2/2]

oomph::GS< CRDoubleMatrix >::GS ( const GS< CRDoubleMatrix > &  )

delete

Broken copy constructor.

Member Function Documentation

clean_up_memory()

void oomph::GS< CRDoubleMatrix >::clean_up_memory ( )

inlineprivatevirtual

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

Reimplemented from oomph::LinearSolver.

    {
      // If the matrix pointer isn't null AND we're allowed to delete the
      // matrix which is only when we create the matrix ourselves
      if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
      {
        // Delete the matrix
        delete Matrix_pt;
 
        // Assign the associated pointer the value NULL
        Matrix_pt = 0;
      }
    } // End of clean_up_memory

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

disable_resolve()

void oomph::GS< CRDoubleMatrix >::disable_resolve ( )

inlinevirtual

Overload disable resolve so that it cleans up memory too.

Reimplemented from oomph::LinearSolver.

    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    } // End of disable_resolve

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

iterations()

unsigned oomph::GS< CRDoubleMatrix >::iterations ( ) const

inlinevirtual

Number of iterations taken.

Implements oomph::IterativeLinearSolver.

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

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

operator=()

void oomph::GS< CRDoubleMatrix >::operator= ( const GS< CRDoubleMatrix > &  )

delete

Broken assignment operator.

preconditioner_setup_time()

double oomph::GS< CRDoubleMatrix >::preconditioner_setup_time ( ) const

inlinevirtual

Returns the time taken to set up the preconditioner.

Reimplemented from oomph::IterativeLinearSolver.

    {
      // Create the error message
      std::string error_output_string = "Gauss Seidel is not a ";
      error_output_string += "preconditionable iterative solver";
 
      // Throw an error
      throw OomphLibError(
        error_output_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
 
      // Return a value so the compiler doesn't throw up an error about
      // no input being returned
      return 0;
    } // End of preconditioner_setup_time

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

resolve()

void oomph::GS< CRDoubleMatrix >::resolve ( const DoubleVector &  rhs, DoubleVector &  result  )

inlinevirtual

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

Reimplemented from oomph::LinearSolver.

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

References oomph::GS< MATRIX >::Matrix_pt, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::GS< MATRIX >::Resolving, and oomph::GS< MATRIX >::solve().

setup_helper()

void oomph::GS< CRDoubleMatrix >::setup_helper

(

DoubleMatrixBase *

matrix_pt

)

Generic setup function to sort out everything that needs to be set up with regards to the input matrix

Explicit template specialisation of the smoother_setup function for CR matrices. Set up the smoother for the matrix specified by the pointer. This definition of the smoother_setup has the added functionality that it sorts the entries in the given matrix so that the CRDoubleMatrix implementation of the solver can be used.

  {
    // Assume the matrix has 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;
 
    // Upcast the input matrix to system matrix to the type MATRIX
    Matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
    // The system matrix here is a CRDoubleMatrix. To make use of the
    // specific implementation of the solver for this type of matrix we
    // need to make sure the entries are arranged correctly
    Matrix_pt->sort_entries();
 
    // Now get access to the vector Index_of_diagonal_entries
    Index_of_diagonal_entries = Matrix_pt->get_index_of_diagonal_entries();
 
#ifdef PARANOID
    // Create a boolean variable to store the result of whether or not
    // the entries are in the correct order
    bool is_matrix_sorted = true;
 
    // Check to make sure the entries have been sorted properly
    is_matrix_sorted = Matrix_pt->entries_are_sorted();
 
    // If the entries are not sorted properly
    if (!is_matrix_sorted)
    {
      // Throw an error; the solver won't work unless the matrix is sorted
      // properly
      throw OomphLibError("Matrix is not sorted correctly",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
  } // End of setup_helper

References oomph::GS< MATRIX >::Matrix_can_be_deleted, oomph::GS< MATRIX >::Matrix_pt, OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

smoother_setup()

void oomph::GS< CRDoubleMatrix >::smoother_setup ( DoubleMatrixBase *  matrix_pt )

inlinevirtual

Set up the smoother for the matrix specified by the pointer.

Implements oomph::Smoother.

    {
      // Assume the matrix has 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;
 
      // Call the generic setup helper function
      setup_helper(matrix_pt);
    } // End of smoother_setup

References oomph::GS< MATRIX >::Matrix_can_be_deleted.

smoother_solve()

void oomph::GS< CRDoubleMatrix >::smoother_solve ( const DoubleVector &  rhs, DoubleVector &  result  )

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

    {
      // If you use a smoother but you don't want to calculate the residual
      Use_as_smoother = true;
 
      // Call the helper function
      solve_helper(Matrix_pt, rhs, result);
    } // End of smoother_solve

References oomph::GS< MATRIX >::Matrix_pt, oomph::GS< MATRIX >::solve_helper(), and oomph::Smoother::Use_as_smoother.

solve() [1/3]

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

    {
      // Reset the Use_as_smoother_flag as the solver is not being used
      // as a smoother
      Use_as_smoother = false;
 
      // Set up the distribution
      this->build_distribution(rhs.distribution_pt());
 
      // Store the matrix if required
      if ((Enable_resolve) && (!Resolving))
      {
        // Upcast to the appropriate matrix type and sort the matrix entries
        // out so that the CRDoubleMatrix implementation of the Gauss-Seidel
        // solver can be used
        Matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
      }
      // We still need to sort the entries
      else
      {
        // The system matrix here is a CRDoubleMatrix. To make use of the
        // specific implementation of the solver for this type of matrix we
        // need to make sure the entries are arranged correctly
        dynamic_cast<CRDoubleMatrix*>(matrix_pt)->sort_entries();
 
        // Now get access to the vector Index_of_diagonal_entries
        Index_of_diagonal_entries = dynamic_cast<CRDoubleMatrix*>(matrix_pt)
                                      ->get_index_of_diagonal_entries();
      }
 
      // Matrix has been passed in from the outside so we must not delete it
      Matrix_can_be_deleted = false;
 
      // Call the helper function
      solve_helper(matrix_pt, rhs, solution);
    } // End of solve

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

solve() [2/3]

void oomph::GS< CRDoubleMatrix >::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() [3/3]

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

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// Explicit template specialisation of the solver for CR matrices: 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.

  {
    // Reset the Use_as_smoother_flag as the solver is not being used
    // as a smoother
    Use_as_smoother = false;
 
    // Find # of degrees of freedom (variables)
    unsigned n_dof = problem_pt->ndof();
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // We're not re-solving
    Resolving = false;
 
    // Get rid of any previously stored data
    clean_up_memory();
 
    // Set up the distribution
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
 
    // Build the distribution of the LinearSolver
    this->build_distribution(dist);
 
    // Allocate space for the Jacobian matrix in CRDoubleMatrix format
    Matrix_pt = new CRDoubleMatrix;
 
    // Build the matrix
    Matrix_pt->build(this->distribution_pt());
 
    // Allocate space for the residuals vector
    DoubleVector f;
 
    // Build it
    f.build(this->distribution_pt(), 0.0);
 
    // Get the global Jacobian and the accompanying residuals vector
    problem_pt->get_jacobian(f, *Matrix_pt);
 
    // Doc. time for setup
    double t_end = TimingHelpers::timer();
    Jacobian_setup_time = t_end - t_start;
 
    // Run the generic setup function
    setup_helper(Matrix_pt);
 
    // We've made the matrix, we can delete it...
    Matrix_can_be_deleted = true;
 
    // If time documentation is enabled
    if (Doc_time)
    {
      // Output the time for the assembly of the Jacobian to screen
      oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
                 << std::endl;
    }
 
    // Call linear algebra-style solver
    solve_helper(Matrix_pt, f, result);
 
    // Kill matrix unless it's still required for resolve
    if (!Enable_resolve) clean_up_memory();
  } // End of solve

References oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::GS< MATRIX >::clean_up_memory(), oomph::Problem::communicator_pt(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Doc_time, oomph::LinearSolver::Enable_resolve, f(), oomph::Problem::get_jacobian(), oomph::IterativeLinearSolver::Jacobian_setup_time, oomph::GS< MATRIX >::Matrix_can_be_deleted, oomph::GS< MATRIX >::Matrix_pt, oomph::Problem::ndof(), oomph::oomph_info, oomph::GS< MATRIX >::Resolving, oomph::GS< MATRIX >::solve_helper(), oomph::TimingHelpers::timer(), and oomph::Smoother::Use_as_smoother.

solve_helper()

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

private

General interface to solve function.

Explicit template specialisation of the solve_helper function for CR matrices. Exploiting the sparsity of the given matrix allows for a much faster iterative solver.

  {
    // Get number of dofs
    unsigned n_dof = rhs.nrow();
 
#ifdef PARANOID
    // Upcast the matrix to the appropriate type
    CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
    // PARANOID Run the self-tests to check the inputs are correct
    this->check_validity_of_solve_helper_inputs<CRDoubleMatrix>(
      cr_matrix_pt, rhs, solution, n_dof);
 
    // We don't need the pointer any more but we do still need the matrix
    // so just make tmp_matrix_pt a null pointer
    cr_matrix_pt = 0;
#endif
 
    // Setup the solution if it is not
    if (!solution.distribution_pt()->built())
    {
      solution.build(this->distribution_pt(), 0.0);
    }
    // If the solution is already set up
    else
    {
      // If we're inside the multigrid solver then as we traverse up the
      // hierarchy we use the smoother on the updated approximate solution.
      // As such, we should ONLY be resetting all the values to zero if
      // we're NOT inside the multigrid solver
      if (!Use_as_smoother)
      {
        // Initialise the vector with all entries set to zero
        solution.initialise(0.0);
      }
    } // if (!solution.distribution_pt()->built())
 
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // Copy the solution vector into x
    DoubleVector x(solution);
 
    // Create a vector to hold the residual. This will only be built if
    // we're not inside the multigrid solver
    DoubleVector current_residual;
 
    // Variable to hold the current residual norm. Only used if we're
    // not inside the multigrid solver
    double norm_res = 0.0;
 
    // Variables to hold the initial residual norm. Only used if we're
    // not inside the multigrid solver
    double norm_f = 0.0;
 
    // Initialise the value of Iterations
    Iterations = 0;
 
    // Calculate the residual only if we're not inside the multigrid solver
    if (!Use_as_smoother)
    {
      // Build the residual vector
      current_residual.build(this->distribution_pt(), 0.0);
 
      // Calculate the residual r=b-Ax
      matrix_pt->residual(x, rhs, current_residual);
 
      // Calculate the 2-norm of the residual vector
      norm_res = current_residual.norm();
 
      // Store the initial norm
      norm_f = norm_res;
 
      // If required will document convergence history to screen
      // or file (if stream is open)
      if (Doc_convergence_history)
      {
        if (!Output_file_stream.is_open())
        {
          oomph_info << Iterations << " " << norm_res << std::endl;
        }
        else
        {
          Output_file_stream << Iterations << " " << norm_res << std::endl;
        }
      } // if (!Use_as_smoother)
    } // if (Doc_convergence_history)
 
    // In each iteration of Gauss-Seidel we need to calculate
    //              x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
    // where L represents the strict lower triangular portion of A,
    // D represents the diagonal of A and U represents the strict upper
    // triangular portion of A. Since the value of inv(L+D)*rhs remains
    // constant throughout we shall calculate it now and store it to
    // avoid calculating it again with each iteration.
 
    // Create a temporary matrix pointer to allow for the extraction of the
    // row start, column index and value pointers
    CRDoubleMatrix* tmp_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
    // First acquire access to the value, row_start and column_index arrays
    // from the compressed row matrix
    const double* value_pt = tmp_matrix_pt->value();
    const int* row_start_pt = tmp_matrix_pt->row_start();
    const int* column_index_pt = tmp_matrix_pt->column_index();
 
    // We've finished using the temporary matrix pointer so make it a null
    // pointer
    tmp_matrix_pt = 0;
 
    // We can only calculate this constant term if the diagonal entries
    // of the matrix are nonzero so we check this first
    for (unsigned i = 0; i < n_dof; i++)
    {
      // Get the index of the last entry below or on the diagonal in row i
      unsigned diag_index = Index_of_diagonal_entries[i];
 
      if (unsigned(*(column_index_pt + diag_index)) != i)
      {
        std::string err_strng = "Gauss-Seidel cannot operate on a matrix with ";
        err_strng += "zero diagonal entries.";
        throw OomphLibError(
          err_strng, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    } // for (unsigned i=0;i<n_dof;i++)
 
    // If there were no zero diagonal entries we can calculate the vector
    //                 constant_term=inv(L+D)*rhs.
 
    // Create a vector to hold the constant term in the iteration
    DoubleVector constant_term(this->distribution_pt(), 0.0);
 
    // Copy the entries of rhs into the vector, constant_term
    constant_term = rhs;
 
    // Start by looping over the rows of constant_term
    for (unsigned i = 0; i < n_dof; i++)
    {
      // Get the index of the last entry below or on the diagonal in row i
      unsigned diag_index = Index_of_diagonal_entries[i];
 
      // Get the index of the first entry in row i
      unsigned row_i_start = *(row_start_pt + i);
 
      // If there are entries strictly below the diagonal then the
      // column index of the first entry in the i-th row cannot be
      // i (since we have ensured that nonzero entries exist on the
      // diagonal so this is the largest the column index of the
      // first entry in this row could be)
      if (unsigned(*(column_index_pt + row_i_start)) != i)
      {
        // Initialise the column index variable
        unsigned column_index = 0;
 
        // Loop over the entries below the diagonal on row i
        for (unsigned j = row_i_start; j < diag_index; j++)
        {
          // Find the column index of this nonzero entry
          column_index = *(column_index_pt + j);
 
          // Subtract the value of A(i,j)*constant_term(j) (where j<i)
          constant_term[i] -= (*(value_pt + j)) * constant_term[column_index];
        }
      } // if (*(column_index_pt+row_i_start)!=i)
 
      // Finish off by calculating constant_term(i)/A(i,i)
      constant_term[i] /= *(value_pt + diag_index);
    } // for (unsigned i=0;i<n_dof;i++)
 
    // Build a temporary vector to store the value of A*x with each iteration
    DoubleVector temp_vec(this->distribution_pt(), 0.0);
 
    // Outermost loop: Run max_iter times (the iteration number)
    for (unsigned iter_num = 0; iter_num < Max_iter; iter_num++)
    {
      // With each iteration we need to calculate
      //            x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
      //                 =-inv(L+D)*U*x_old+constant_term.
      // Since we've already calculated the vector constant_term, it remains
      // for us to calculate inv(L+D)*U*x. To do this we break the calculation
      // into two parts:
      //        (1) The matrix-vector multiplication temp_vec=U*x, and;
      //        (2) The matrix-vector multiplication inv(L+D)*temp_vec.
 
      //------------------------------
      // (1) Calculating temp_vec=U*x:
      //------------------------------
 
      // Auxiliary variable to store the index of the first entry on the
      // upper triangular portion of the matrix
      unsigned upper_tri_start = 0;
 
      // Auxiliary variable to store the index of the first entry in the
      // next row
      unsigned next_row_start = 0;
 
      // Set the temporary vector temp_vec to be initially be the zero vector
      temp_vec.initialise(0.0);
 
      // Loop over the rows of temp_vec (note: we do not loop over the final
      // row since all the entries on the last row of a strict upper triangular
      // matrix are zero)
      for (unsigned i = 0; i < n_dof - 1; i++)
      {
        // Get the index of the first entry on the upper triangular portion of
        // the matrix noting that the i-th entry in Index_of_diagonal_entries
        // will store the index of the diagonal entry of the i-th row
        upper_tri_start = Index_of_diagonal_entries[i] + 1;
 
        // Get the index of the first entry in the next row
        next_row_start = *(row_start_pt + i + 1);
 
        // Variable to store the column index of each entry
        unsigned column_index = 0;
 
        // Loop over all of the entries above the diagonal
        for (unsigned j = upper_tri_start; j < next_row_start; j++)
        {
          // Get the column index of this entry
          column_index = *(column_index_pt + j);
 
          // Update the value of temp_vec[i] by adding A(i,j)*x(j) (where j>i)
          temp_vec[i] += (*(value_pt + j)) * x[column_index];
        }
      } // for (unsigned i=0;i<n_dof;i++)
 
      //-----------------------------------
      // (2) Calculating inv(L+D)*temp_vec:
      //-----------------------------------
 
      // Now U*x=temp_vec has been calculated we need to calculate the vector
      //                 y=inv(L+D)*U*x=inv(L+D)*temp_vec.
      // Note: the i-th entry in y is given by
      //       y(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*y(j))/A(i,i).
      // We can see from this that to calculate the current entry in y we
      // use its previously calculated entries and the current entry in
      // temp_vec. As a result, we do not need two separate vectors y and
      // temp_vec for this calculation. Instead, we can just update the
      // entries in temp_vec (carefully!) using:
      //   temp_vec(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*temp_vec(j))/A(i,i).
 
      // Start by looping over the rows of rhs
      for (unsigned i = 0; i < n_dof; i++)
      {
        // Get the index of the last entry below or on the diagonal in row i
        unsigned diag_index = Index_of_diagonal_entries[i];
 
        // Get the index of the first entry in row i
        unsigned row_i_start = *(row_start_pt + i);
 
        // If there are no entries strictly below the diagonal then the
        // column index of the first entry in the i-th row will be greater
        // than or equal to i
        if (unsigned(*(column_index_pt + row_i_start)) != i)
        {
          // Initialise the column index variable
          unsigned column_index = 0;
 
          // Loop over the entries below the diagonal
          for (unsigned j = row_i_start; j < diag_index; j++)
          {
            // Find the column index of this nonzero entry
            column_index = *(column_index_pt + j);
 
            // Subtract the value of A(i,j)*y(j) (where j<i)
            temp_vec[i] -= (*(value_pt + j)) * temp_vec[column_index];
          }
        } // if (*(column_index_pt+row_i_start)!=i)
 
        // Finish off by calculating temp_vec(i)/A(i,i)
        temp_vec[i] /= *(value_pt + diag_index);
      } // for (unsigned i=0;i<n_dof;i++)
 
      //-------------------------
      // Updating the solution x:
      //-------------------------
 
      // The updated solution is given by
      //            x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
      //                 =-temp_vec+constant_term,
      // We have calculated both constant_term and temp_vec so we simply
      // need to calculate x now
 
      // Set x to be the vector, constant_term
      x = constant_term;
 
      // Update x to give the next iterate
      x -= temp_vec;
 
      // Increment the value of Iterations
      Iterations++;
 
      // Calculate the residual only if we're not inside the multigrid solver
      if (!Use_as_smoother)
      {
        // Get residual
        matrix_pt->residual(x, rhs, current_residual);
 
        // Calculate the relative residual norm (i.e.
        // \frac{\|r_{i}\|}{\|r_{0}\|})
        norm_res = current_residual.norm() / norm_f;
 
        // 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_num << " " << norm_res << std::endl;
          }
          else
          {
            Output_file_stream << iter_num << " " << norm_res << std::endl;
          }
        } // if (Doc_convergence_history)
 
        // Check the tolerance only if the residual norm is being computed
        if (norm_res < Tolerance)
        {
          // Break out of the for-loop
          break;
        }
      } // if (!Use_as_smoother)
    } // for (unsigned iter_num=0;iter_num<max_iter;iter_num++)
 
    // Calculate the residual only if we're not inside the multigrid solver
    if (!Use_as_smoother)
    {
      // If time documentation is enabled
      if (Doc_time)
      {
        oomph_info << "\nGS converged. Residual norm: " << norm_res
                   << "\nNumber of iterations to convergence: " << Iterations
                   << "\n"
                   << std::endl;
      }
    } // if (!Use_as_smoother)
 
    // Copy the solution into the solution vector
    solution = x;
 
    // Document the time for the solver
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
 
    // If time documentation is enabled
    if (Doc_time)
    {
      oomph_info << "Time for solve with Gauss-Seidel [sec]: " << Solution_time
                 << std::endl;
    }
 
    // If the solver failed to converge and the user asked for an error if
    // this happened
    if ((Iterations > Max_iter - 1) && (Throw_error_after_max_iter))
    {
      std::string error_message =
        "Solver failed to converge and you requested ";
      error_message += "an error on convergence failures.";
      throw OomphLibError(
        error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
  } // End of solve_helper

References oomph::DoubleVector::build(), oomph::CRDoubleMatrix::column_index(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::IterativeLinearSolver::Doc_convergence_history, oomph::LinearSolver::Doc_time, i, oomph::DoubleVector::initialise(), oomph::GS< MATRIX >::Iterations, j, oomph::IterativeLinearSolver::Max_iter, oomph::DoubleVector::norm(), oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::IterativeLinearSolver::Output_file_stream, oomph::DoubleMatrixBase::residual(), oomph::CRDoubleMatrix::row_start(), BiharmonicTestFunctions1::solution(), oomph::IterativeLinearSolver::Solution_time, oomph::Global_string_for_annotation::string(), oomph::IterativeLinearSolver::Throw_error_after_max_iter, oomph::TimingHelpers::timer(), oomph::IterativeLinearSolver::Tolerance, oomph::Smoother::Use_as_smoother, oomph::CRDoubleMatrix::value(), and plotDoE::x.

Member Data Documentation

Index_of_diagonal_entries

Vector<int> oomph::GS< CRDoubleMatrix >::Index_of_diagonal_entries

private

Vector whose i'th entry contains the index of the last entry below or on the diagonal of the i'th row of the matrix

Iterations

unsigned oomph::GS< CRDoubleMatrix >::Iterations

private

Number of iterations taken.

Matrix_can_be_deleted

bool oomph::GS< CRDoubleMatrix >::Matrix_can_be_deleted

private

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

Matrix_pt

CRDoubleMatrix* oomph::GS< CRDoubleMatrix >::Matrix_pt

private

System matrix pointer in the format specified by the template argument.

Resolving

bool oomph::GS< CRDoubleMatrix >::Resolving

private

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

---

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

• iterative_linear_solver.h
• iterative_linear_solver.cc