#include <matrices.h>

Inheritance diagram for oomph::CRDoubleMatrix:

Classes
struct	CRDoubleMatrixComparisonHelper
	Create a struct to provide a comparison function for std::sort. More...

Public Member Functions
	CRDoubleMatrix ()
	Default constructor. More...

	CRDoubleMatrix (const LinearAlgebraDistribution *distribution_pt, const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)

	CRDoubleMatrix (const LinearAlgebraDistribution *distribution_pt)

	CRDoubleMatrix (const CRDoubleMatrix &matrix)
	Copy constructor. More...

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

virtual	~CRDoubleMatrix ()
	Destructor. More...

const Vector< int >	get_index_of_diagonal_entries () const

bool	entries_are_sorted (const bool &doc_unordered_entries=false) const

void	sort_entries ()

void	build (const LinearAlgebraDistribution *distribution_pt, const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)

void	build (const LinearAlgebraDistribution *distribution_pt)
	Rebuild the matrix - assembles an empty matrix with a defined distribution. More...

void	build (const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)
	keeps the existing distribution and just matrix that is stored More...

void	build_without_copy (const unsigned &ncol, const unsigned &nnz, double value, int column_index, int *row_start)
	method to rebuild the matrix, but not the distribution More...

void	redistribute (const LinearAlgebraDistribution *const &dist_pt)

void	clear ()
	clear More...

unsigned long	nrow () const
	Return the number of rows of the matrix. More...

unsigned long	ncol () const
	Return the number of columns of the matrix. More...

void	output_bottom_right_zero_helper (std::ostream &outfile) const

void	sparse_indexed_output_helper (std::ostream &outfile) const

void	sparse_indexed_output_with_offset (std::string filename)

double	operator() (const unsigned long &i, const unsigned long &j) const

int *	row_start ()
	Access to C-style row_start array. More...

const int *	row_start () const
	Access to C-style row_start array (const version) More...

int *	column_index ()
	Access to C-style column index array. More...

const int *	column_index () const
	Access to C-style column index array (const version) More...

double *	value ()
	Access to C-style value array. More...

const double *	value () const
	Access to C-style value array (const version) More...

unsigned long	nnz () const
	Return the number of nonzero entries (the local nnz) More...

virtual void	ludecompose ()
	Do LU decomposition. More...

virtual void	lubksub (DoubleVector &rhs)
	LU back solve for given RHS. More...

void	multiply (const DoubleVector &x, DoubleVector &soln) const
	Multiply the matrix by the vector x: soln=Ax. More...

void	multiply_transpose (const DoubleVector &x, DoubleVector &soln) const
	Multiply the transposed matrix by the vector x: soln=A^T x. More...

void	multiply (const CRDoubleMatrix &matrix_in, CRDoubleMatrix &result) const

void	matrix_reduction (const double &alpha, CRDoubleMatrix &reduced_matrix)

unsigned &	serial_matrix_matrix_multiply_method ()

const unsigned &	serial_matrix_matrix_multiply_method () const

unsigned &	distributed_matrix_matrix_multiply_method ()

const unsigned &	distributed_matrix_matrix_multiply_method () const

bool	built () const

CRDoubleMatrix *	global_matrix () const

void	get_matrix_transpose (CRDoubleMatrix *result) const
	Returns the transpose of this matrix. More...

double	inf_norm () const
	returns the inf-norm of this matrix More...

Vector< double >	diagonal_entries () const

void	add (const CRDoubleMatrix &matrix_in, CRDoubleMatrix &result_matrix) const
	element-wise addition of this matrix with matrix_in. More...

Public Member Functions inherited from oomph::Matrix< double, CRDoubleMatrix >
	Matrix ()
	(Empty) constructor More...

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

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

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

double	operator() (const unsigned long &i, const unsigned long &j) const

double &	operator() (const unsigned long &i, const unsigned long &j)

virtual void	output (std::ostream &outfile) const

void	sparse_indexed_output (std::ostream &outfile, const unsigned &precision=0, const bool &output_bottom_right_zero=false) const

void	sparse_indexed_output (std::string filename, const unsigned &precision=0, const bool &output_bottom_right_zero=false) const

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

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

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

virtual	~DoubleMatrixBase ()
	virtual (empty) destructor More...

LinearSolver *&	linear_solver_pt ()
	Return a pointer to the linear solver object. More...

LinearSolver *const &	linear_solver_pt () const
	Return a pointer to the linear solver object (const version) More...

void	solve (DoubleVector &rhs)

void	solve (const DoubleVector &rhs, DoubleVector &soln)

void	solve (Vector< double > &rhs)

void	solve (const Vector< double > &rhs, Vector< double > &soln)

virtual void	residual (const DoubleVector &x, const DoubleVector &b, DoubleVector &residual_)
	Find the residual, i.e. r=b-Ax the residual. More...

virtual double	max_residual (const DoubleVector &x, const DoubleVector &rhs)

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

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

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

virtual	~DistributableLinearAlgebraObject ()
	Destructor. More...

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

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

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

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

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

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

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

bool	distribution_built () const

void	build_distribution (const LinearAlgebraDistribution *const dist_pt)

void	build_distribution (const LinearAlgebraDistribution &dist)

Public Attributes
struct oomph::CRDoubleMatrix::CRDoubleMatrixComparisonHelper	Comparison_struct

Private Attributes
Vector< int >	Index_of_diagonal_entries

unsigned	Serial_matrix_matrix_multiply_method

unsigned	Distributed_matrix_matrix_multiply_method

CRMatrix< double >	CR_matrix
	Storage for the Matrix in CR Format. More...

bool	Built

Additional Inherited Members
Protected Member Functions inherited from oomph::Matrix< double, CRDoubleMatrix >
void	range_check (const unsigned long &i, const unsigned long &j) const

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

Protected Attributes inherited from oomph::DoubleMatrixBase
LinearSolver *	Linear_solver_pt

LinearSolver *	Default_linear_solver_pt

Detailed Description

A class for compressed row matrices. This is a distributable object.

Constructor & Destructor Documentation

◆ CRDoubleMatrix() [1/4]

oomph::CRDoubleMatrix::CRDoubleMatrix ( )

Default constructor.

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

   {
     // set the default solver
     Linear_solver_pt = Default_linear_solver_pt = new SuperLUSolver;
  
     // matrix not built
     Built = false;
  
     // set the serial matrix-matrix multiply method
 #ifdef OOMPH_HAS_TRILINOS
     //    Serial_matrix_matrix_multiply_method = 4;
     Serial_matrix_matrix_multiply_method = 2;
 #else
     Serial_matrix_matrix_multiply_method = 2;
 #endif
   }

References Built, oomph::DoubleMatrixBase::Default_linear_solver_pt, oomph::DoubleMatrixBase::Linear_solver_pt, and Serial_matrix_matrix_multiply_method.

Referenced by global_matrix().

◆ CRDoubleMatrix() [2/4]

oomph::CRDoubleMatrix::CRDoubleMatrix	(	const LinearAlgebraDistribution *	dist_pt,
		const unsigned &	ncol,
		const Vector< double > &	value,
		const Vector< int > &	column_index,
		const Vector< int > &	row_start
	)

Constructor: vector of values, vector of column indices, vector of row starts and number of rows and columns.

Constructor: Takes the distribution and the number of columns, as well as the vector of values, vector of column indices,vector of row starts.

   {
     // build the compressed row matrix
     CR_matrix.build(
       value, column_index, row_start, dist_pt->nrow_local(), ncol);
  
     // store the Distribution
     this->build_distribution(dist_pt);
  
     // set the linear solver
     Linear_solver_pt = Default_linear_solver_pt = new SuperLUSolver;
  
     // set the serial matrix-matrix multiply method
 #ifdef OOMPH_HAS_TRILINOS
     // Serial_matrix_matrix_multiply_method = 4;
     Serial_matrix_matrix_multiply_method = 2;
 #else
     Serial_matrix_matrix_multiply_method = 2;
 #endif
  
     // matrix has been built
     Built = true;
   }

References oomph::CRMatrix< T >::build(), oomph::DistributableLinearAlgebraObject::build_distribution(), Built, column_index(), CR_matrix, oomph::DoubleMatrixBase::Default_linear_solver_pt, oomph::DoubleMatrixBase::Linear_solver_pt, ncol(), oomph::LinearAlgebraDistribution::nrow_local(), row_start(), Serial_matrix_matrix_multiply_method, and value().

◆ CRDoubleMatrix() [3/4]

oomph::CRDoubleMatrix::CRDoubleMatrix ( const LinearAlgebraDistribution * distribution_pt )

Constructor: just stores the distribution but does not build the matrix

   {
     this->build_distribution(distribution_pt);
  
     // set the default solver
     Linear_solver_pt = Default_linear_solver_pt = new SuperLUSolver;
  
     // matrix not built
     Built = false;
  
 // set the serial matrix-matrix multiply method
 #ifdef OOMPH_HAS_TRILINOS
     //    Serial_matrix_matrix_multiply_method = 4;
     Serial_matrix_matrix_multiply_method = 2;
 #else
     Serial_matrix_matrix_multiply_method = 2;
 #endif
   }

References oomph::DistributableLinearAlgebraObject::build_distribution(), Built, oomph::DoubleMatrixBase::Default_linear_solver_pt, oomph::DoubleMatrixBase::Linear_solver_pt, and Serial_matrix_matrix_multiply_method.

◆ CRDoubleMatrix() [4/4]

oomph::CRDoubleMatrix::CRDoubleMatrix ( const CRDoubleMatrix & matrix )

Copy constructor.

   {
     // copy the distribution
     this->build_distribution(other_matrix.distribution_pt());
  
     // copy coefficients
     const double* values_pt = other_matrix.value();
     const int* column_indices = other_matrix.column_index();
     const int* row_start = other_matrix.row_start();
  
     // This is the local nnz.
     const unsigned nnz = other_matrix.nnz();
  
     // Using number of local rows since the underlying CRMatrix is local to
     // each processor.
     const unsigned nrow_local = other_matrix.nrow_local();
  
     // Storage for the (yet to be copied) data.
     double* my_values_pt = new double[nnz];
     int* my_column_indices = new int[nnz];
     int* my_row_start = new int[nrow_local + 1];
  
     // Copying over the data.
     std::copy(values_pt, values_pt + nnz, my_values_pt);
     std::copy(column_indices, column_indices + nnz, my_column_indices);
     std::copy(row_start, row_start + nrow_local + 1, my_row_start);
  
  
     // Build without copy since we have made a deep copy of the data structure.
     this->build_without_copy(
       other_matrix.ncol(), nnz, my_values_pt, my_column_indices, my_row_start);
  
     // set the default solver
     Linear_solver_pt = Default_linear_solver_pt = new SuperLUSolver;
  
     // matrix is built
     Built = true;
  
     // set the serial matrix-matrix multiply method
 #ifdef OOMPH_HAS_TRILINOS
     // Serial_matrix_matrix_multiply_method = 4;
     Serial_matrix_matrix_multiply_method = 2;
 #else
     Serial_matrix_matrix_multiply_method = 2;
 #endif
   }

References oomph::DistributableLinearAlgebraObject::build_distribution(), build_without_copy(), Built, column_index(), copy(), oomph::DoubleMatrixBase::Default_linear_solver_pt, oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::DoubleMatrixBase::Linear_solver_pt, ncol(), nnz(), oomph::DistributableLinearAlgebraObject::nrow_local(), row_start(), Serial_matrix_matrix_multiply_method, and value().

◆ ~CRDoubleMatrix()

oomph::CRDoubleMatrix::~CRDoubleMatrix ( )

virtual

Destructor.

   {
     this->clear();
     delete Default_linear_solver_pt;
     Default_linear_solver_pt = 0;
   }

References clear(), and oomph::DoubleMatrixBase::Default_linear_solver_pt.

Member Function Documentation

◆ add()

void oomph::CRDoubleMatrix::add	(	const CRDoubleMatrix &	matrix_in,
		CRDoubleMatrix &	result_matrix
	)		const

element-wise addition of this matrix with matrix_in.

Element-wise addition of this matrix with matrix_in.

   {
 #ifdef PARANOID
     // Check if this matrix is built.
     if (!this->built())
     {
       std::ostringstream error_message;
       error_message << "The matrix is not built.\n"
                     << "Please build the matrix!\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // Check if this matrix_in is built.
     if (!matrix_in.built())
     {
       std::ostringstream error_message;
       error_message << "The matrix matrix_in is not built.\n"
                     << "Please build the matrix!\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // Check if the dimensions of this matrix and matrix_in are the same.
     unsigned long this_nrow = this->nrow();
     unsigned long matrix_in_nrow = matrix_in.nrow();
     if (this_nrow != matrix_in_nrow)
     {
       std::ostringstream error_message;
       error_message << "matrix_in has a different number of rows than\n"
                     << "this matrix.\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     unsigned long this_ncol = this->ncol();
     unsigned long matrix_in_ncol = matrix_in.ncol();
     if (this_ncol != matrix_in_ncol)
     {
       std::ostringstream error_message;
       error_message << "matrix_in has a different number of columns than\n"
                     << "this matrix.\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // Check if the distribution is the same (Otherwise we may have to send and
     // receive data from other processors - which is not implemented!)
     if (*(this->distribution_pt()) != *(matrix_in.distribution_pt()))
     {
       std::ostringstream error_message;
       error_message << "matrix_in must have the same distribution as\n"
                     << "this matrix.\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // If the matrix is built, check that it's existing distribution is the
     // same as the in matrix. Since we'll use the same distribution instead
     // of completely rebuilding it.
     if (result_matrix.built() &&
         (*result_matrix.distribution_pt() != *matrix_in.distribution_pt()))
     {
       std::ostringstream error_message;
       error_message << "The result_matrix is built. "
                     << "But has a different distribution from matrix_in \n"
                     << "They need to be the same.\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // To add the elements of two CRDoubleMatrices, we need to know the union of
     // the sparsity patterns. This is determined by the column indices.
     // We add the column indices and values (entries) as a key-value pair in
     // a map (per row). We then read these out into two column indices and
     // values vector for the result matrix.
  
     unsigned nrow_local = this->nrow_local();
     Vector<int> res_column_indices;
     Vector<double> res_values;
     Vector<int> res_row_start;
     res_row_start.reserve(nrow_local + 1);
  
     // The row_start and column_indices
     const int* this_column_indices = this->column_index();
     const int* this_row_start = this->row_start();
     const int* in_column_indices = matrix_in.column_index();
     const int* in_row_start = matrix_in.row_start();
  
     // Values from this matrix and matrix_in.
     const double* this_values = this->value();
     const double* in_values = matrix_in.value();
  
  
     // The first entry in row_start is always zero.
     res_row_start.push_back(0);
  
     // Loop through the rows of both matrices and insert the column indices and
     // values as a key-value pair.
     for (unsigned row_i = 0; row_i < nrow_local; row_i++)
     {
       // Create the map for this row.
       std::map<int, double> res_row_map;
  
       // Insert the column and value pair for this matrix.
       for (int i = this_row_start[row_i]; i < this_row_start[row_i + 1]; i++)
       {
         res_row_map[this_column_indices[i]] = this_values[i];
       }
  
       // Insert the column and value pair for in matrix.
       for (int i = in_row_start[row_i]; i < in_row_start[row_i + 1]; i++)
       {
         res_row_map[in_column_indices[i]] += in_values[i];
       }
  
       // Fill in the row start
       res_row_start.push_back(res_row_start.back() + res_row_map.size());
  
       // Now insert the key into res_column_indices and value into res_values
       for (std::map<int, double>::iterator it = res_row_map.begin();
            it != res_row_map.end();
            ++it)
       {
         res_column_indices.push_back(it->first);
         res_values.push_back(it->second);
       }
     }
  
     // Finally build the result_matrix.
     if (result_matrix.distribution_pt()->built())
     {
       // Use the existing distribution.
       result_matrix.build(
         this->ncol(), res_values, res_column_indices, res_row_start);
     }
     else
     {
       // Build with THIS distribution
       result_matrix.build(this->distribution_pt(),
                           this->ncol(),
                           res_values,
                           res_column_indices,
                           res_row_start);
     }
   }

References build(), oomph::LinearAlgebraDistribution::built(), built(), column_index(), oomph::DistributableLinearAlgebraObject::distribution_pt(), i, ncol(), nrow(), oomph::DistributableLinearAlgebraObject::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, row_start(), and value().

Referenced by main().

◆ build() [1/3]

void oomph::CRDoubleMatrix::build ( const LinearAlgebraDistribution * distribution_pt )

Rebuild the matrix - assembles an empty matrix with a defined distribution.

rebuild the matrix - assembles an empty matrix will a defined distribution

   {
     this->clear();
     this->build_distribution(distribution_pt);
   }

References oomph::DistributableLinearAlgebraObject::build_distribution(), and clear().

◆ build() [2/3]

void oomph::CRDoubleMatrix::build	(	const LinearAlgebraDistribution *	distribution_pt,
		const unsigned &	ncol,
		const Vector< double > &	value,
		const Vector< int > &	column_index,
		const Vector< int > &	row_start
	)

build method: vector of values, vector of column indices, vector of row starts and number of rows and columns.

build method: Takes the distribution and the number of columns, as well as the vector of values, vector of column indices,vector of row starts.

   {
     // clear
     this->clear();
  
     // store the Distribution
     this->build_distribution(distribution_pt);
  
     // set the linear solver
     Default_linear_solver_pt = new SuperLUSolver;
  
     // now build the matrix
     this->build(ncol, value, column_index, row_start);
   }

References oomph::DistributableLinearAlgebraObject::build_distribution(), clear(), column_index(), oomph::DoubleMatrixBase::Default_linear_solver_pt, row_start(), and value().

Referenced by add(), oomph::CRDoubleMatrixHelpers::concatenate(), oomph::CRDoubleMatrixHelpers::concatenate_without_communication(), oomph::CRDoubleMatrixHelpers::create_uniformly_distributed_matrix(), oomph::CRDoubleMatrixHelpers::deep_copy(), oomph::Problem::get_eigenproblem_matrices(), oomph::Problem::get_jacobian(), get_matrix_transpose(), oomph::BlockPreconditioner< MATRIX >::internal_get_block(), main(), matrix_reduction(), oomph::TrilinosEpetraHelpers::multiply(), multiply(), redistribute(), oomph::LagrangeEnforcedFlowPreconditioner::setup(), and oomph::HelmholtzMGPreconditioner< DIM >::setup_mg_structures().

◆ build() [3/3]

void oomph::CRDoubleMatrix::build	(	const unsigned &	ncol,
		const Vector< double > &	value,
		const Vector< int > &	column_index,
		const Vector< int > &	row_start
	)

keeps the existing distribution and just matrix that is stored

method to rebuild the matrix, but not the distribution

   {
     // call the underlying build method
     CR_matrix.clean_up_memory();
     CR_matrix.build(value, column_index, row_start, this->nrow_local(), ncol);
  
     // matrix has been build
     Built = true;
   }

References oomph::CRMatrix< T >::build(), Built, oomph::CRMatrix< T >::clean_up_memory(), column_index(), CR_matrix, oomph::DistributableLinearAlgebraObject::nrow_local(), row_start(), and value().

◆ build_without_copy()

void oomph::CRDoubleMatrix::build_without_copy	(	const unsigned &	ncol,
		const unsigned &	nnz,
		double *	value,
		int *	column_index,
		int *	row_start
	)

method to rebuild the matrix, but not the distribution

keeps the existing distribution and just matrix that is stored without copying the matrix data

   {
     // call the underlying build method
     CR_matrix.clean_up_memory();
     CR_matrix.build_without_copy(
       value, column_index, row_start, nnz, this->nrow_local(), ncol);
  
     // matrix has been build
     Built = true;
   }

References oomph::CRMatrix< T >::build_without_copy(), Built, oomph::CRMatrix< T >::clean_up_memory(), column_index(), CR_matrix, nnz(), oomph::DistributableLinearAlgebraObject::nrow_local(), row_start(), and value().

Referenced by oomph::NavierStokesSchurComplementPreconditioner::assemble_inv_press_and_veloc_mass_matrix_diagonal(), oomph::PressureBasedSolidLSCPreconditioner::assemble_mass_matrix_diagonal(), oomph::CRDoubleMatrixHelpers::concatenate_without_communication(), CRDoubleMatrix(), oomph::CRDoubleMatrixHelpers::deep_copy(), oomph::Problem::get_eigenproblem_matrices(), oomph::Problem::get_jacobian(), global_matrix(), oomph::BlockPreconditioner< MATRIX >::internal_get_block(), oomph::TrilinosEpetraHelpers::multiply(), multiply(), and redistribute().

◆ built()

bool oomph::CRDoubleMatrix::built ( ) const

inline

access function to the Built flag - indicates whether the matrix has been build - i.e. the distribution has been defined and the matrix assembled.

     {
       return Built;
     }

References Built.

Referenced by add(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix_for_aztecoo(), oomph::HypreHelpers::create_HYPRE_Matrix(), oomph::CRDoubleMatrixHelpers::create_uniformly_distributed_matrix(), oomph::CRDoubleMatrixHelpers::deep_copy(), diagonal_entries(), oomph::MumpsSolver::factorise(), oomph::BlockPreconditioner< MATRIX >::get_concatenated_block(), multiply(), oomph::TrilinosEpetraHelpers::multiply(), oomph::DoubleVector::norm(), and redistribute().

◆ clear()

void oomph::CRDoubleMatrix::clear ( )

clear

Clean method.

   {
     this->clear_distribution();
     CR_matrix.clean_up_memory();
     Built = false;
  
     if (Linear_solver_pt != 0) // Only clean up if it exists
       Linear_solver_pt->clean_up_memory();
   }

References Built, oomph::LinearSolver::clean_up_memory(), oomph::CRMatrix< T >::clean_up_memory(), oomph::DistributableLinearAlgebraObject::clear_distribution(), CR_matrix, and oomph::DoubleMatrixBase::Linear_solver_pt.

Referenced by build(), oomph::MumpsSolver::factorise(), oomph::HypreInterface::hypre_matrix_setup(), main(), oomph::FSIPreconditioner::setup(), and ~CRDoubleMatrix().

◆ column_index() [1/2]

int* oomph::CRDoubleMatrix::column_index ( )

inline

Access to C-style column index array.

     {
       return CR_matrix.column_index();
     }

References oomph::CRMatrix< T >::column_index(), and CR_matrix.

Referenced by add(), build(), build_without_copy(), oomph::CRDoubleMatrixHelpers::concatenate(), CRDoubleMatrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix_for_aztecoo(), oomph::HypreHelpers::create_HYPRE_Matrix(), oomph::CRDoubleMatrixHelpers::deep_copy(), entries_are_sorted(), oomph::MumpsSolver::factorise(), oomph::SuperLUSolver::factorise_serial(), get_matrix_transpose(), global_matrix(), oomph::BlockPreconditioner< MATRIX >::internal_get_block(), matrix_reduction(), multiply(), multiply_transpose(), redistribute(), oomph::ILUZeroPreconditioner< CRDoubleMatrix >::setup(), oomph::PseudoElasticPreconditionerSubsidiaryPreconditionerOld::setup(), oomph::HelmholtzMGPreconditioner< DIM >::setup_coarsest_level_structures(), oomph::GS< CRDoubleMatrix >::solve_helper(), sort_entries(), and sparse_indexed_output_with_offset().

◆ column_index() [2/2]

const int* oomph::CRDoubleMatrix::column_index ( ) const

inline

Access to C-style column index array (const version)

     {
       return CR_matrix.column_index();
     }

References oomph::CRMatrix< T >::column_index(), and CR_matrix.

◆ diagonal_entries()

Vector< double > oomph::CRDoubleMatrix::diagonal_entries ( ) const

returns a Vector of diagonal entries of this matrix. This only works with square matrices. This condition may be relaxed in the future if need be.

Return the diagonal entries of the matrix. This only works with square matrices. This condition may be relaxed in the future if need be.

   {
 #ifdef PARANOID
     // Check if the matrix has been built.
     if (!this->built())
     {
       std::ostringstream error_message;
       error_message << "The matrix has not been built.\n"
                     << "Please build it...\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // Check if this is a square matrix.
     if (this->nrow() != this->ncol())
     {
       std::ostringstream error_message;
       error_message << "The matrix is not square. Can only get the diagonal\n"
                     << "entries of a square matrix.\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // We get the diagonal entries on this processor.
     unsigned nrow_local = this->nrow_local();
  
     // Create storage for the diagonal entries.
     Vector<double> result_vec;
     result_vec.reserve(nrow_local);
  
     // Get the first row for the column offset.
     unsigned first_row = this->first_row();
  
     // Loop through the local rows.
     for (unsigned i = 0; i < nrow_local; i++)
     {
       // The column entries are globally indexed.
       unsigned diag_entry_col = first_row + i;
  
       // Push back the diagonal entry.
       result_vec.push_back(CR_matrix.get_entry(i, diag_entry_col));
     }
  
     return result_vec;
   }

References built(), CR_matrix, oomph::DistributableLinearAlgebraObject::first_row(), oomph::CRMatrix< T >::get_entry(), i, ncol(), nrow(), oomph::DistributableLinearAlgebraObject::nrow_local(), OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

Referenced by oomph::ComplexDampedJacobi< MATRIX >::complex_smoother_setup(), and oomph::LagrangeEnforcedFlowPreconditioner::setup().

◆ distributed_matrix_matrix_multiply_method() [1/2]

unsigned& oomph::CRDoubleMatrix::distributed_matrix_matrix_multiply_method ( )

inline

Access function to Distributed_matrix_matrix_multiply_method, the flag which determines the matrix matrix multiplication method used for distributed matrices. Method 1: Trilinos Epetra Matrix Matrix multiply. Method 2: Trilinos Epetra Matrix Matrix multiply (ML based).

     {
       return Distributed_matrix_matrix_multiply_method;
     }

References Distributed_matrix_matrix_multiply_method.

Referenced by oomph::CRDoubleMatrixHelpers::deep_copy().

◆ distributed_matrix_matrix_multiply_method() [2/2]

const unsigned& oomph::CRDoubleMatrix::distributed_matrix_matrix_multiply_method ( ) const

inline

Read only access function (const version) to Distributed_matrix_matrix_multiply_method, the flag which determines the matrix matrix multiplication method used for distributed matrices. Method 1: Trilinos Epetra Matrix Matrix multiply. Method 2: Trilinos Epetra Matrix Matrix multiply (ML based).

     {
       return Distributed_matrix_matrix_multiply_method;
     }

References Distributed_matrix_matrix_multiply_method.

◆ entries_are_sorted()

bool oomph::CRDoubleMatrix::entries_are_sorted ( const bool & doc_unordered_entries = false ) const

Runs through the column index vector and checks if the entries follow the regular lexicographical ordering of matrix entries, i.e. it will check (at the i-th row of the matrix) if the entries in the column index vector associated with this row are in increasing order

Runs through the column index vector and checks if the entries are arranged arbitrarily or if they follow the regular lexicographical of matrices. If a boolean argument is provided with the assignment TRUE then information on the first entry which is not in the correct position will also be given

   {
 #ifdef OOMPH_HAS_MPI
     // We only need to produce a warning if the matrix is distributed
     if (this->distributed())
     {
       // Create an ostringstream object to store the warning message
       std::ostringstream warning_message;
  
       // Create the warning messsage
       warning_message << "This method currently works for serial but "
                       << "has not been implemented for use with MPI!\n";
  
       // Issue the warning
       OomphLibWarning(warning_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Get the number of rows in this matrix
     unsigned n_rows = this->nrow();
  
     // Acquire access to the value, row_start and column_index arrays from
     // the CR matrix. Since we do not change anything in row_start_pt we
     // give it the const prefix
     const int* column_index_pt = this->column_index();
     const int* row_start_pt = this->row_start();
  
     // Loop over the rows of matrix
     for (unsigned i = 0; i < n_rows; i++)
     {
       // Calculate the number of nonzeros in the i-th row
       unsigned nnz_row_i = *(row_start_pt + i + 1) - *(row_start_pt + i);
  
       // Get the index of the first entry in row i
       unsigned i_row_start = *(row_start_pt + i);
  
       // Loop over the entries of the i-th row
       for (unsigned j = 0; j < nnz_row_i - 1; j++)
       {
         // Check if the column index of the following entry is greater than the
         // current entry
         if ((*(column_index_pt + i_row_start + j + 1)) <
             (*(column_index_pt + i_row_start + j)))
         {
           // If the input argument was set to TRUE we document we output
           // information of the first entry which is not in the correct position
           if (doc_unordered_entries)
           {
             // Tell the user
             oomph_info << "Matrix has not been correctly sorted!"
                        << "\nOn row: " << i << "\nEntry: " << j
                        << "\nEntry 1 = " << *(column_index_pt + i_row_start + j)
                        << "\nEntry 2 = "
                        << *(column_index_pt + i_row_start + j + 1) << std::endl;
           }
  
           // It hasn't worked
           return false;
         } // if ((*(column_index_pt+i_row_start+j+1)) ...
       } // for (unsigned j=0;j<nnz_row_i-1;j++)
     } // for (unsigned i=0;i<n_rows;i++)
  
     // If it gets here without a warning then the entries in each row of
     // the matrix are ordered by increasing column index
     return true;
   } // End of entries_are_sorted()

References column_index(), oomph::DistributableLinearAlgebraObject::distributed(), i, j, nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, and row_start().

◆ get_index_of_diagonal_entries()

const Vector<int> oomph::CRDoubleMatrix::get_index_of_diagonal_entries ( ) const

inline

Access function: returns the vector Index_of_diagonal_entries. The i-th entry of the vector contains the index of the last entry below or on the diagonal. If there are no entries below or on the diagonal then the corresponding entry is -1. If, however, there are no entries in the row then the entry is irrelevant and is kept as the initialised value; 0.

     {
       // Check to see if the vector has been set up
       if (Index_of_diagonal_entries.size() == 0)
       {
         // Make the warning
         std::string err_strng =
           "The Index_of_diagonal_entries vector has not been ";
         err_strng += "set up yet. Run sort_entries() to set this vector up.";
  
         // Throw the warning
         OomphLibWarning(err_strng,
                         "CRDoubleMatrix::get_index_of_diagonal_entries()",
                         OOMPH_EXCEPTION_LOCATION);
       }
  
       // Return the vector
       return Index_of_diagonal_entries;
     } // End of index_of_diagonal_entries

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

◆ get_matrix_transpose()

void oomph::CRDoubleMatrix::get_matrix_transpose ( CRDoubleMatrix * result ) const

Returns the transpose of this matrix.

Compute transpose of matrix.

   {
     // Get the number of non_zeros
     unsigned long nnon_zeros = this->nnz();
  
     // Find the number of rows and columns in the transposed
     // matrix. We differentiate these from those associated
     // with the original matrix by appending the characters
     // '_t' onto the end
     unsigned long n_rows = this->nrow();
     unsigned long n_rows_t = this->ncol();
  
 #ifdef OOMPH_HAS_MPI
     // We only need to produce a warning if the matrix is distributed
     if (this->distributed())
     {
       // Create an ostringstream object to store the warning message
       std::ostringstream warning_message;
  
       // Create the warning messsage
       warning_message << "This method currently works for serial but "
                       << "has not been tested with MPI!\n";
  
       // Issue the warning
       OomphLibWarning(warning_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Set up the distribution for the transposed matrix
     result->distribution_pt()->build(
       this->distribution_pt()->communicator_pt(), n_rows_t, false);
  
     // Acquire access to the value, row_start and column_index
     // arrays from the CR matrix
     const double* value_pt = this->value();
     const int* row_start_pt = this->row_start();
     const int* column_index_pt = this->column_index();
  
     // Allocate space for the row_start and column_index vectors
     // associated with the transpose of the current matrix.
     Vector<double> value_t(nnon_zeros, 0.0);
     Vector<int> column_index_t(nnon_zeros, 0);
     Vector<int> row_start_t(n_rows_t + 1, 0);
  
     // Loop over the column index vector and count how many times
     // each column number occurs and increment the appropriate
     // entry in row_start_t whose i+1'th entry will contain the
     // number of non-zeros in the i'th column of the matrix (this
     // is done so that the cumulative sum done next returns the
     // correct result)
     for (unsigned i = 0; i < nnon_zeros; i++)
     {
       // Assign entries to row_start_t (noting the first
       // entry is left as 0 for the cumulative sum done later)
       row_start_t[*(column_index_pt + i) + 1]++;
     }
  
     // Calculate the sum of the first i entries in the row_start_t
     // vector and store the values in the i'th entry of row_start_t
     for (unsigned i = 1; i < n_rows_t + 1; i++)
     {
       // Calculate the cumulative sum
       row_start_t[i] += row_start_t[i - 1];
     }
  
     // Allocate space for variables to store the indices of the
     // start and end of the non-zeros in a given row of the matrix
     unsigned i_row_s = 0;
     unsigned i_row_e = 0;
  
     // Initialise 3 extra variables for readability of the
     // code in the subsequent piece of code
     unsigned a = 0;
     unsigned b = 0;
     unsigned c = 0;
  
     // Vector needed to count the number of entries added
     // to each segment in column_index_t where each segment
     // is associated with each row in the transposed matrix
     Vector<int> counter(n_rows_t, 0);
  
     // Set the entries in column_index_t. To do this we loop
     // over each row of the original matrix (equivalently
     // the number of columns in the transpose)
     for (unsigned i_row = 0; i_row < n_rows; i_row++)
     {
       // Here we find the indices of the start and end
       // of the non-zeros in i_row'th row of the matrix.
       // [Note, there should be a -1 on i_row_e but this
       // is ignored so that we use a strict inequality
       // in the subsequent for-loop!]
       i_row_s = *(row_start_pt + i_row);
       i_row_e = *(row_start_pt + i_row + 1);
  
       // Loop over the entries in the i_row'th row
       // of the matrix
       for (unsigned j = i_row_s; j < i_row_e; j++)
       {
         // The value of a is the column index of the j'th
         // element in the i_row'th row of the original matrix
         // (which is also the row index of the associated
         // non-zero in the transposed matrix)
         a = *(column_index_pt + j);
  
         // The value of b will be used to find the start
         // of the appropriate segment of column_index_t
         // that the non-zero belongs in
         b = row_start_t[a];
  
         // Find the number of elements already added to
         // this segment (to find which entry of the segment
         // the value i_row, the column index of the non-zero
         // in the transposed matrix, should be assigned to)
         c = counter[*(column_index_pt + j)];
  
         // Assign the value i_row to the appropriate entry
         // in column_index_t
         column_index_t[b + c] = i_row;
         value_t[b + c] = *(value_pt + j);
  
         // Increment the j'th value of the vector counter
         // to indicate another non-zero index has been
         // added into the
         counter[*(column_index_pt + j)]++;
  
       } // End of for-loop over i_row'th row of the matrix
  
     } // End of for-loop over rows of the matrix
  
     // Build the matrix (note: the value of n_cols for the
     // transposed matrix is n_rows for the original matrix)
     result->build(n_rows, value_t, column_index_t, row_start_t);
  
   } // End of the function

References a, b, build(), oomph::LinearAlgebraDistribution::build(), calibrate::c, column_index(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), i, j, ncol(), nnz(), nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, row_start(), and value().

◆ global_matrix()

CRDoubleMatrix * oomph::CRDoubleMatrix::global_matrix ( ) const

if this matrix is distributed then a the equivalent global matrix is built using new and returned. The calling method is responsible for the destruction of the new matrix.

if this matrix is distributed then the equivalent global matrix is built using new and returned. The calling method is responsible for the destruction of the new matrix.

   {
 #ifdef OOMPH_HAS_MPI
     // if this matrix is not distributed then this method is redundant
     if (!this->distributed() ||
         this->distribution_pt()->communicator_pt()->nproc() == 1)
     {
       return new CRDoubleMatrix(*this);
     }
  
     // nnz
     int nnz = this->nnz();
  
     // my nrow local
     unsigned nrow_local = this->nrow_local();
  
     // nrow global
     unsigned nrow = this->nrow();
  
     // cache nproc
     int nproc = this->distribution_pt()->communicator_pt()->nproc();
  
     // get the nnzs on the other processors
     int* dist_nnz_pt = new int[nproc];
     MPI_Allgather(&nnz,
                   1,
                   MPI_INT,
                   dist_nnz_pt,
                   1,
                   MPI_INT,
                   this->distribution_pt()->communicator_pt()->mpi_comm());
  
     // create a int vector of first rows and nrow local and compute nnz global
     int* dist_first_row = new int[nproc];
     int* dist_nrow_local = new int[nproc];
     int nnz_global = 0;
     for (int p = 0; p < nproc; p++)
     {
       nnz_global += dist_nnz_pt[p];
       dist_first_row[p] = this->first_row(p);
       dist_nrow_local[p] = this->nrow_local(p);
     }
  
     // compute the offset for the values and column index data
     int* nnz_offset = new int[nproc];
     nnz_offset[0] = 0;
     for (int p = 1; p < nproc; p++)
     {
       nnz_offset[p] = nnz_offset[p - 1] + dist_nnz_pt[p - 1];
     }
  
     // get pointers to the (current) distributed data
     // const_cast required because MPI requires non-const data when sending
     // data
     int* dist_row_start = const_cast<int*>(this->row_start());
     int* dist_column_index = const_cast<int*>(this->column_index());
     double* dist_value = const_cast<double*>(this->value());
  
     // space for the global matrix
     int* global_row_start = new int[nrow + 1];
     int* global_column_index = new int[nnz_global];
     double* global_value = new double[nnz_global];
  
     // get the row starts
     MPI_Allgatherv(dist_row_start,
                    nrow_local,
                    MPI_INT,
                    global_row_start,
                    dist_nrow_local,
                    dist_first_row,
                    MPI_INT,
                    this->distribution_pt()->communicator_pt()->mpi_comm());
  
     // get the column indexes
     MPI_Allgatherv(dist_column_index,
                    nnz,
                    MPI_INT,
                    global_column_index,
                    dist_nnz_pt,
                    nnz_offset,
                    MPI_INT,
                    this->distribution_pt()->communicator_pt()->mpi_comm());
  
     // get the values
     MPI_Allgatherv(dist_value,
                    nnz,
                    MPI_DOUBLE,
                    global_value,
                    dist_nnz_pt,
                    nnz_offset,
                    MPI_DOUBLE,
                    this->distribution_pt()->communicator_pt()->mpi_comm());
  
     // finally the last row start
     global_row_start[nrow] = nnz_global;
  
     // update the other row start
     for (int p = 0; p < nproc; p++)
     {
       for (int i = 0; i < dist_nrow_local[p]; i++)
       {
         unsigned j = dist_first_row[p] + i;
         global_row_start[j] += nnz_offset[p];
       }
     }
  
     // create the global distribution
     LinearAlgebraDistribution* dist_pt = new LinearAlgebraDistribution(
       this->distribution_pt()->communicator_pt(), nrow, false);
  
     // create the matrix
     CRDoubleMatrix* matrix_pt = new CRDoubleMatrix(dist_pt);
  
     // copy of distribution taken so delete
     delete dist_pt;
  
     // pass data into matrix
     matrix_pt->build_without_copy(this->ncol(),
                                   nnz_global,
                                   global_value,
                                   global_column_index,
                                   global_row_start);
  
     // clean up
     delete[] dist_first_row;
     delete[] dist_nrow_local;
     delete[] nnz_offset;
     delete[] dist_nnz_pt;
  
     // and return
     return matrix_pt;
 #else
     return new CRDoubleMatrix(*this);
 #endif
   }

References build_without_copy(), column_index(), oomph::LinearAlgebraDistribution::communicator_pt(), CRDoubleMatrix(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::DistributableLinearAlgebraObject::first_row(), i, j, ncol(), nnz(), oomph::OomphCommunicator::nproc(), nrow(), oomph::DistributableLinearAlgebraObject::nrow_local(), p, row_start(), and value().

Referenced by oomph::ILUZeroPreconditioner< CRDoubleMatrix >::setup().

◆ inf_norm()

double oomph::CRDoubleMatrix::inf_norm ( ) const

returns the inf-norm of this matrix

Compute infinity (maximum) norm of matrix.

   {
 #ifdef PARANOID
     // paranoid check that the vector is setup
     if (!this->distribution_built())
     {
       std::ostringstream error_message;
       error_message << "This matrix must be setup.";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // compute the local norm
     unsigned nrow_local = this->nrow_local();
     double n = 0;
     const int* row_start = CR_matrix.row_start();
     const double* value = CR_matrix.value();
     for (unsigned i = 0; i < nrow_local; i++)
     {
       double a = 0;
       for (int j = row_start[i]; j < row_start[i + 1]; j++)
       {
         a += fabs(value[j]);
       }
       n = std::max(n, a);
     }
  
     // if this vector is distributed and on multiple processors then gather
 #ifdef OOMPH_HAS_MPI
     double n2 = n;
     if (this->distributed() &&
         this->distribution_pt()->communicator_pt()->nproc() > 1)
     {
       MPI_Allreduce(&n,
                     &n2,
                     1,
                     MPI_DOUBLE,
                     MPI_MAX,
                     this->distribution_pt()->communicator_pt()->mpi_comm());
     }
     n = n2;
 #endif
  
     // and return
     return n;
   }

References a, CR_matrix, oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_built(), oomph::DistributableLinearAlgebraObject::distribution_pt(), boost::multiprecision::fabs(), i, j, max, n, oomph::DistributableLinearAlgebraObject::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRMatrix< T >::row_start(), row_start(), oomph::SparseMatrix< T, MATRIX_TYPE >::value(), and value().

Referenced by oomph::PseudoElasticPreconditionerScalingHelperOld::s_inf_norm().

◆ lubksub()

void oomph::CRDoubleMatrix::lubksub ( DoubleVector & rhs )

virtual

LU back solve for given RHS.

Do back-substitution.

   {
 #ifdef PARANOID
     // check that the rhs vector is setup
     if (!rhs.built())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The vector rhs has not been setup";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // check that the rhs vector has the same distribution as this matrix
     if (!(*this->distribution_pt() == *rhs.distribution_pt()))
     {
       std::ostringstream error_message_stream;
       error_message_stream
         << "The vector rhs must have the same distribution as the matrix";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // backsub
     DoubleVector rhs_copy(rhs);
     static_cast<SuperLUSolver*>(Default_linear_solver_pt)
       ->backsub(rhs_copy, rhs);
   }

References oomph::DoubleVector::built(), oomph::DoubleMatrixBase::Default_linear_solver_pt, oomph::DistributableLinearAlgebraObject::distribution_pt(), OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ ludecompose()

void oomph::CRDoubleMatrix::ludecompose ( )

virtual

Do LU decomposition.

LU decomposition using SuperLU if matrix is not distributed or distributed onto a single processor.

   {
 #ifdef PARANOID
     // check that the this matrix is built
     if (!Built)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "This matrix has not been built.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // factorise using superlu or superlu dist if we oomph has mpi
     static_cast<SuperLUSolver*>(Default_linear_solver_pt)->factorise(this);
   }

References Built, oomph::DoubleMatrixBase::Default_linear_solver_pt, OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ matrix_reduction()

void oomph::CRDoubleMatrix::matrix_reduction	(	const double &	alpha,
		CRDoubleMatrix &	reduced_matrix
	)

For every row, find the maximum absolute value of the entries in this row. Set all values that are less than alpha times this maximum to zero and return the resulting matrix in reduced_matrix. Note: Diagonal entries are retained regardless of their size.

   {
     // number of rows in matrix
     long n_row = nrow_local();
     double max_row;
  
     // Here's the packed format for the new matrix
     Vector<int> B_row_start(1);
     Vector<int> B_column_index;
     Vector<double> B_value;
  
     // get pointers to the underlying data
     const int* row_start = CR_matrix.row_start();
     const int* column_index = CR_matrix.column_index();
     const double* value = CR_matrix.value();
  
     // k is counter for the number of entries in the reduced matrix
     unsigned k = 0;
  
     // Initialise row start
     B_row_start[0] = 0;
  
     // Loop over rows
     for (long i = 0; i < n_row; i++)
     {
       // Initialise max value in row
       max_row = 0.0;
  
       // Loop over entries in columns
       for (long j = row_start[i]; j < row_start[i + 1]; j++)
       {
         // Find max. value in row
         if (std::fabs(value[j]) > max_row)
         {
           max_row = std::fabs(value[j]);
         }
       }
  
       // Decide if we need to retain the entries in the row
       for (long j = row_start[i]; j < row_start[i + 1]; j++)
       {
         // If we're on the diagonal or the value is sufficiently large: retain
         // i.e. copy across.
         if (i == column_index[j] || std::fabs(value[j]) > alpha * max_row)
         {
           B_value.push_back(value[j]);
           B_column_index.push_back(column_index[j]);
           k++;
         }
       }
       // This writes the row start for the next row -- equal to
       // to the number of entries written so far (C++ zero-based indexing!)
       B_row_start.push_back(k);
     }
  
     // Build the matrix from the compressed format
     dynamic_cast<CRDoubleMatrix&>(reduced_matrix)
       .build(this->ncol(), B_value, B_column_index, B_row_start);
   }

References alpha, build(), oomph::CRMatrix< T >::column_index(), column_index(), CR_matrix, boost::multiprecision::fabs(), i, j, k, ncol(), oomph::DistributableLinearAlgebraObject::nrow_local(), oomph::CRMatrix< T >::row_start(), row_start(), oomph::SparseMatrix< T, MATRIX_TYPE >::value(), and value().

◆ multiply() [1/2]

void oomph::CRDoubleMatrix::multiply	(	const CRDoubleMatrix &	matrix_in,
		CRDoubleMatrix &	result
	)		const

Function to multiply this matrix by the CRDoubleMatrix matrix_in. In a serial matrix, there are 4 methods available: Method 1: First runs through this matrix and matrix_in to find the storage requirements for result - arrays of the correct size are then allocated before performing the calculation. Minimises memory requirements but more costly. Method 2: Grows storage for values and column indices of result 'on the fly' using an array of maps. Faster but more memory intensive. Method 3: Grows storage for values and column indices of result 'on the fly' using a vector of vectors. Not particularly impressive on the platforms we tried... Method 4: Trilinos Epetra Matrix Matrix multiply. Method 5: Trilinox Epetra Matrix Matrix Mulitply (ml based) If Trilinos is installed then Method 4 is employed by default, otherwise Method 2 is employed by default. In a distributed matrix, only Trilinos Epetra Matrix Matrix multiply is available.

   {
 #ifdef PARANOID
     // check that this matrix is built
     if (!Built)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "This matrix has not been built";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // check that this matrix is built
     if (!matrix_in.built())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "This matrix matrix_in has not been built";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // if soln is setup then it should have the same distribution as x
     if (result.built())
     {
       if (!(*result.distribution_pt() == *this->distribution_pt()))
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "The matrix result is setup and therefore must have the same "
           << "distribution as the vector x";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // if the result has not been setup, then store the distribution
     if (!result.distribution_built())
     {
       result.build(this->distribution_pt());
     }
  
     // short name for Serial_matrix_matrix_multiply_method
     unsigned method = Serial_matrix_matrix_multiply_method;
  
     // if this matrix is not distributed and matrix in is not distributed
     if (!this->distributed() && !matrix_in.distributed() &&
         ((method == 1) || (method == 2) || (method == 3)))
     {
       // NB N is number of rows!
       unsigned long N = this->nrow();
       unsigned long M = matrix_in.ncol();
       unsigned long Nnz = 0;
  
       // pointers to arrays which store result
       int* Row_start = 0;
       double* Value = 0;
       int* Column_index = 0;
  
       // get pointers to matrix_in
       const int* matrix_in_row_start = matrix_in.row_start();
       const int* matrix_in_column_index = matrix_in.column_index();
       const double* matrix_in_value = matrix_in.value();
  
       // get pointers to this matrix
       const double* this_value = this->value();
       const int* this_row_start = this->row_start();
       const int* this_column_index = this->column_index();
  
       // clock_t clock1 = clock();
  
       // METHOD 1
       // --------
       if (method == 1)
       {
         // allocate storage for row starts
         Row_start = new int[N + 1];
         Row_start[0] = 0;
  
         // a set to store number of non-zero columns in each row of result
         std::set<unsigned> columns;
  
         // run through rows of this matrix and matrix_in to find number of
         // non-zero entries in each row of result
         for (unsigned long this_row = 0; this_row < N; this_row++)
         {
           // run through non-zeros in this_row of this matrix
           for (int this_ptr = this_row_start[this_row];
                this_ptr < this_row_start[this_row + 1];
                this_ptr++)
           {
             // find column index for non-zero
             int matrix_in_row = this_column_index[this_ptr];
  
             // run through corresponding row in matrix_in
             for (int matrix_in_ptr = matrix_in_row_start[matrix_in_row];
                  matrix_in_ptr < matrix_in_row_start[matrix_in_row + 1];
                  matrix_in_ptr++)
             {
               // find column index for non-zero in matrix_in and store in
               // columns
               columns.insert(matrix_in_column_index[matrix_in_ptr]);
             }
           }
           // update Row_start
           Row_start[this_row + 1] = Row_start[this_row] + columns.size();
  
           // wipe values in columns
           columns.clear();
         }
  
         // set Nnz
         Nnz = Row_start[N];
  
         // allocate arrays for result
         Value = new double[Nnz];
         Column_index = new int[Nnz];
  
         // set all values of Column_index to -1
         for (unsigned long i = 0; i < Nnz; i++)
         {
           Column_index[i] = -1;
         }
  
         // Calculate values for result - first run through rows of this matrix
         for (unsigned long this_row = 0; this_row < N; this_row++)
         {
           // run through non-zeros in this_row
           for (int this_ptr = this_row_start[this_row];
                this_ptr < this_row_start[this_row + 1];
                this_ptr++)
           {
             // find value of non-zero
             double this_val = this_value[this_ptr];
  
             // find column associated with non-zero
             int matrix_in_row = this_column_index[this_ptr];
  
             // run through corresponding row in matrix_in
             for (int matrix_in_ptr = matrix_in_row_start[matrix_in_row];
                  matrix_in_ptr < matrix_in_row_start[matrix_in_row + 1];
                  matrix_in_ptr++)
             {
               // find column index for non-zero in matrix_in
               int col = matrix_in_column_index[matrix_in_ptr];
  
               // find position in result to insert value
               for (int ptr = Row_start[this_row];
                    ptr <= Row_start[this_row + 1];
                    ptr++)
               {
                 if (ptr == Row_start[this_row + 1])
                 {
                   // error - have passed end of row without finding
                   // correct column
                   std::ostringstream error_message;
                   error_message << "Error inserting value in result";
  
                   throw OomphLibError(error_message.str(),
                                       OOMPH_CURRENT_FUNCTION,
                                       OOMPH_EXCEPTION_LOCATION);
                 }
                 else if (Column_index[ptr] == -1)
                 {
                   // first entry for this column index
                   Column_index[ptr] = col;
                   Value[ptr] = this_val * matrix_in_value[matrix_in_ptr];
                   break;
                 }
                 else if (Column_index[ptr] == col)
                 {
                   // column index already exists - add value
                   Value[ptr] += this_val * matrix_in_value[matrix_in_ptr];
                   break;
                 }
               }
             }
           }
         }
       }
  
       // METHOD 2
       // --------
       else if (method == 2)
       {
         // generate array of maps to store values for result
         std::map<int, double>* result_maps = new std::map<int, double>[N];
  
         // run through rows of this matrix
         for (unsigned long this_row = 0; this_row < N; this_row++)
         {
           // run through non-zeros in this_row
           for (int this_ptr = this_row_start[this_row];
                this_ptr < this_row_start[this_row + 1];
                this_ptr++)
           {
             // find value of non-zero
             double this_val = this_value[this_ptr];
  
             // find column index associated with non-zero
             int matrix_in_row = this_column_index[this_ptr];
  
             // run through corresponding row in matrix_in
             for (int matrix_in_ptr = matrix_in_row_start[matrix_in_row];
                  matrix_in_ptr < matrix_in_row_start[matrix_in_row + 1];
                  matrix_in_ptr++)
             {
               // find column index for non-zero in matrix_in
               int col = matrix_in_column_index[matrix_in_ptr];
  
               // insert value
               result_maps[this_row][col] +=
                 this_val * matrix_in_value[matrix_in_ptr];
             }
           }
         }
  
         // allocate Row_start
         Row_start = new int[N + 1];
  
         // copy across row starts
         Row_start[0] = 0;
         for (unsigned long row = 0; row < N; row++)
         {
           int size = result_maps[row].size();
           Row_start[row + 1] = Row_start[row] + size;
         }
  
         // set Nnz
         Nnz = Row_start[N];
  
         // allocate other arrays
         Value = new double[Nnz];
         Column_index = new int[Nnz];
  
         // copy values and column indices
         for (unsigned long row = 0; row < N; row++)
         {
           unsigned ptr = Row_start[row];
           for (std::map<int, double>::iterator i = result_maps[row].begin();
                i != result_maps[row].end();
                i++)
           {
             Column_index[ptr] = i->first;
             Value[ptr] = i->second;
             ptr++;
           }
         }
  
         // tidy up memory
         delete[] result_maps;
       }
  
       // METHOD 3
       // --------
       else if (method == 3)
       {
         // vectors of vectors to store results
         std::vector<std::vector<int>> result_cols(N);
         std::vector<std::vector<double>> result_vals(N);
  
         // run through the rows of this matrix
         for (unsigned long this_row = 0; this_row < N; this_row++)
         {
           // run through non-zeros in this_row
           for (int this_ptr = this_row_start[this_row];
                this_ptr < this_row_start[this_row + 1];
                this_ptr++)
           {
             // find value of non-zero
             double this_val = this_value[this_ptr];
  
             // find column index associated with non-zero
             int matrix_in_row = this_column_index[this_ptr];
  
             // run through corresponding row in matrix_in
             for (int matrix_in_ptr = matrix_in_row_start[matrix_in_row];
                  matrix_in_ptr < matrix_in_row_start[matrix_in_row + 1];
                  matrix_in_ptr++)
             {
               // find column index for non-zero in matrix_in
               int col = matrix_in_column_index[matrix_in_ptr];
  
               // insert value
               int size = result_cols[this_row].size();
               for (int i = 0; i <= size; i++)
               {
                 if (i == size)
                 {
                   // first entry for this column
                   result_cols[this_row].push_back(col);
                   result_vals[this_row].push_back(
                     this_val * matrix_in_value[matrix_in_ptr]);
                 }
                 else if (col == result_cols[this_row][i])
                 {
                   // column already exists
                   result_vals[this_row][i] +=
                     this_val * matrix_in_value[matrix_in_ptr];
                   break;
                 }
               }
             }
           }
         }
  
         // allocate Row_start
         Row_start = new int[N + 1];
  
         // copy across row starts
         Row_start[0] = 0;
         for (unsigned long row = 0; row < N; row++)
         {
           int size = result_cols[row].size();
           Row_start[row + 1] = Row_start[row] + size;
         }
  
         // set Nnz
         Nnz = Row_start[N];
  
         // allocate other arrays
         Value = new double[Nnz];
         Column_index = new int[Nnz];
  
         // copy across values and column indices
         for (unsigned long row = 0; row < N; row++)
         {
           unsigned ptr = Row_start[row];
           unsigned nnn = result_cols[row].size();
           for (unsigned i = 0; i < nnn; i++)
           {
             Column_index[ptr] = result_cols[row][i];
             Value[ptr] = result_vals[row][i];
             ptr++;
           }
         }
       }
  
       // build
       result.build_without_copy(M, Nnz, Value, Column_index, Row_start);
     }
  
     // else we have to use trilinos
     else
     {
 #ifdef OOMPH_HAS_TRILINOS
       bool use_ml = false;
       if (method == 5)
       {
         use_ml = true;
       }
       TrilinosEpetraHelpers::multiply(*this, matrix_in, result, use_ml);
 #else
       std::ostringstream error_message;
       error_message << "Serial_matrix_matrix_multiply_method = "
                     << Serial_matrix_matrix_multiply_method
                     << " requires trilinos.";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
 #endif
     }
   }

References build(), build_without_copy(), built(), Built, col(), column_index(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_built(), oomph::DistributableLinearAlgebraObject::distribution_pt(), i, oomph::TrilinosEpetraHelpers::multiply(), N, ncol(), nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, row(), row_start(), Serial_matrix_matrix_multiply_method, size, and value().

◆ multiply() [2/2]

void oomph::CRDoubleMatrix::multiply	(	const DoubleVector &	x,
		DoubleVector &	soln
	)		const

virtual

Multiply the matrix by the vector x: soln=Ax.

Multiply the matrix by the vector x.

Implements oomph::DoubleMatrixBase.

   {
 #ifdef PARANOID
     // check that this matrix is built
     if (!Built)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "This matrix has not been built";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // check that the distribution of x is setup
     if (!x.built())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The distribution of the vector x must be setup";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // Check to see if x.size() = ncol().
     if (this->ncol() != x.distribution_pt()->nrow())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The number of rows in the x vector and the "
                               "number of columns in the "
                            << "matrix must be the same";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // if the soln is distributed
     if (soln.built())
     {
       if (!(*soln.distribution_pt() == *this->distribution_pt()))
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "The soln vector is setup and therefore must have the same "
           << "distribution as the matrix";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // if soln is not setup then setup the distribution
     if (!soln.built())
     {
       // Resize and initialize the solution vector
       soln.build(this->distribution_pt(), 0.0);
     }
  
     // Initialise
     soln.initialise(0.0);
  
     // if distributed and on more than one processor use trilinos
     // otherwise use the oomph-lib methods
     if (this->distributed() &&
         this->distribution_pt()->communicator_pt()->nproc() > 1)
     {
 #ifdef OOMPH_HAS_TRILINOS
       // This will only work if we have trilinos on board
       TrilinosEpetraHelpers::multiply(this, x, soln);
 #else
       std::ostringstream error_message_stream;
       error_message_stream
         << "Matrix-vector product on multiple processors with distributed "
         << "CRDoubleMatrix requires Trilinos.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
 #endif
     }
     else
     {
       unsigned n = this->nrow();
       const int* row_start = CR_matrix.row_start();
       const int* column_index = CR_matrix.column_index();
       const double* value = CR_matrix.value();
       double* soln_pt = soln.values_pt();
       const double* x_pt = x.values_pt();
       for (unsigned long i = 0; i < n; i++)
       {
         soln_pt[i] = 0.0;
         for (long k = row_start[i]; k < row_start[i + 1]; k++)
         {
           unsigned long j = column_index[k];
           double a_ij = value[k];
           soln_pt[i] += a_ij * x_pt[j];
         }
       }
     }
   }

References oomph::DoubleVector::build(), oomph::DoubleVector::built(), Built, oomph::CRMatrix< T >::column_index(), column_index(), CR_matrix, oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), i, oomph::DoubleVector::initialise(), j, k, oomph::TrilinosEpetraHelpers::multiply(), n, ncol(), oomph::LinearAlgebraDistribution::nrow(), nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRMatrix< T >::row_start(), row_start(), oomph::SparseMatrix< T, MATRIX_TYPE >::value(), value(), oomph::DoubleVector::values_pt(), and plotDoE::x.

Referenced by oomph::ProblemBasedShiftInvertOperator::apply(), oomph::AugmentedProblemGMRES::augmented_matrix_multiply(), main(), oomph::MatrixVectorProduct::multiply(), oomph::DoubleVector::norm(), oomph::GMRESBlockPreconditioner::preconditioner_solve(), oomph::LagrangeEnforcedFlowPreconditioner::setup(), oomph::NavierStokesSchurComplementPreconditioner::setup(), oomph::PressureBasedSolidLSCPreconditioner::setup(), and oomph::SpaceTimeNavierStokesSubsidiaryPreconditioner::setup().

◆ multiply_transpose()

void oomph::CRDoubleMatrix::multiply_transpose	(	const DoubleVector &	x,
		DoubleVector &	soln
	)		const

virtual

Multiply the transposed matrix by the vector x: soln=A^T x.

Implements oomph::DoubleMatrixBase.

   {
 #ifdef PARANOID
     // check that this matrix is built
     if (!Built)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "This matrix has not been built";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // Check to see if x.size() = ncol().
     if (!(*this->distribution_pt() == *x.distribution_pt()))
     {
       std::ostringstream error_message_stream;
       error_message_stream
         << "The x vector and this matrix must have the same distribution.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // if soln is setup then it should have the same distribution as x
     if (soln.built())
     {
       if (soln.distribution_pt()->nrow() != this->ncol())
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "The soln vector is setup and therefore must have the same "
           << "number of rows as the vector x";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // if soln is not setup then setup the distribution
     if (!soln.built())
     {
       LinearAlgebraDistribution* dist_pt =
         new LinearAlgebraDistribution(x.distribution_pt()->communicator_pt(),
                                       this->ncol(),
                                       this->distributed());
       soln.build(dist_pt, 0.0);
       delete dist_pt;
     }
  
     // Initialise
     soln.initialise(0.0);
  
     if (this->distributed() &&
         this->distribution_pt()->communicator_pt()->nproc() > 1)
     {
 #ifdef OOMPH_HAS_TRILINOS
       // This will only work if we have trilinos on board
       TrilinosEpetraHelpers::multiply(this, x, soln);
 #else
       std::ostringstream error_message_stream;
       error_message_stream
         << "Matrix-vector product on multiple processors with distributed "
         << "CRDoubleMatrix requires Trilinos.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
 #endif
     }
     else
     {
       unsigned n = this->nrow();
       const int* row_start = CR_matrix.row_start();
       const int* column_index = CR_matrix.column_index();
       const double* value = CR_matrix.value();
       double* soln_pt = soln.values_pt();
       const double* x_pt = x.values_pt();
       // Matrix vector product
       for (unsigned long i = 0; i < n; i++)
       {
         for (long k = row_start[i]; k < row_start[i + 1]; k++)
         {
           unsigned long j = column_index[k];
           double a_ij = value[k];
           soln_pt[j] += a_ij * x_pt[i];
         }
       }
     }
   }

References oomph::DoubleVector::build(), oomph::DoubleVector::built(), Built, oomph::CRMatrix< T >::column_index(), column_index(), CR_matrix, oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), i, oomph::DoubleVector::initialise(), j, k, oomph::TrilinosEpetraHelpers::multiply(), n, oomph::LinearAlgebraDistribution::nrow(), nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRMatrix< T >::row_start(), row_start(), oomph::SparseMatrix< T, MATRIX_TYPE >::value(), value(), oomph::DoubleVector::values_pt(), and plotDoE::x.

Referenced by oomph::MatrixVectorProduct::multiply_transpose(), oomph::SuperLUSolver::solve(), and oomph::SuperLUSolver::solve_transpose().

◆ ncol()

unsigned long oomph::CRDoubleMatrix::ncol ( ) const

inlinevirtual

Return the number of columns of the matrix.

Implements oomph::DoubleMatrixBase.

     {
       return CR_matrix.ncol();
     }

References CR_matrix, and oomph::SparseMatrix< T, MATRIX_TYPE >::ncol().

Referenced by add(), CRDoubleMatrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix_for_aztecoo(), oomph::HypreHelpers::create_HYPRE_Matrix(), oomph::CRDoubleMatrixHelpers::deep_copy(), diagonal_entries(), oomph::BlockPreconditioner< MATRIX >::get_concatenated_block(), get_matrix_transpose(), global_matrix(), main(), matrix_reduction(), oomph::TrilinosEpetraHelpers::multiply(), multiply(), redistribute(), oomph::MatrixVectorProduct::setup(), and oomph::HelmholtzMGPreconditioner< DIM >::setup_coarsest_level_structures().

◆ nnz()

unsigned long oomph::CRDoubleMatrix::nnz ( ) const

inline

Return the number of nonzero entries (the local nnz)

     {
       return CR_matrix.nnz();
     }

References CR_matrix, and oomph::SparseMatrix< T, MATRIX_TYPE >::nnz().

Referenced by build_without_copy(), CRDoubleMatrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix_for_aztecoo(), oomph::CRDoubleMatrixHelpers::deep_copy(), oomph::MumpsSolver::factorise(), oomph::SuperLUSolver::factorise_serial(), get_matrix_transpose(), global_matrix(), main(), redistribute(), oomph::HyprePreconditioner::setup(), oomph::LagrangeEnforcedFlowPreconditioner::setup(), oomph::BandedBlockTriangularPreconditioner< MATRIX >::setup(), oomph::SpaceTimeNavierStokesSubsidiaryPreconditioner::setup(), oomph::GMRESBlockPreconditioner::setup(), oomph::HelmholtzMGPreconditioner< DIM >::setup_coarsest_level_structures(), oomph::SuperLUSolver::solve(), and oomph::SuperLUSolver::solve_transpose().

◆ nrow()

unsigned long oomph::CRDoubleMatrix::nrow ( ) const

inlinevirtual

Return the number of rows of the matrix.

Implements oomph::DoubleMatrixBase.

     {
       return DistributableLinearAlgebraObject::nrow();
     }

References oomph::DistributableLinearAlgebraObject::nrow().

Referenced by add(), oomph::AugmentedProblemGMRES::augmented_matrix_multiply(), oomph::ComplexDampedJacobi< MATRIX >::complex_smoother_setup(), oomph::CRDoubleMatrixHelpers::concatenate(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix_for_aztecoo(), oomph::HypreHelpers::create_HYPRE_Matrix(), diagonal_entries(), entries_are_sorted(), oomph::Problem::get_eigenproblem_matrices(), get_matrix_transpose(), global_matrix(), oomph::HypreInterface::hypre_matrix_setup(), main(), oomph::MGSolver< DIM >::modify_restriction_matrices(), oomph::TrilinosEpetraHelpers::multiply(), multiply(), multiply_transpose(), redistribute(), oomph::ILUZeroPreconditioner< CRDoubleMatrix >::setup(), oomph::HyprePreconditioner::setup(), oomph::BandedBlockTriangularPreconditioner< MATRIX >::setup(), oomph::SpaceTimeNavierStokesSubsidiaryPreconditioner::setup(), oomph::GMRESBlockPreconditioner::setup(), oomph::HelmholtzMGPreconditioner< DIM >::setup_coarsest_level_structures(), oomph::SuperLUSolver::solve(), oomph::SuperLUSolver::solve_transpose(), and sort_entries().

◆ operator()()

double oomph::CRDoubleMatrix::operator()	(	const unsigned long &	i,
		const unsigned long &	j
	)		const

inlinevirtual

Overload the round-bracket access operator for read-only access. In a distributed matrix i refers to the local row index.

Implements oomph::DoubleMatrixBase.

     {
       return CR_matrix.get_entry(i, j);
     }

References CR_matrix, oomph::CRMatrix< T >::get_entry(), i, and j.

◆ operator=()

void oomph::CRDoubleMatrix::operator= ( const CRDoubleMatrix & )

delete

Broken assignment operator.

◆ output_bottom_right_zero_helper()

void oomph::CRDoubleMatrix::output_bottom_right_zero_helper ( std::ostream & outfile ) const

inlinevirtual

Output the "bottom right" entry regardless of it being zero or not (this allows automatic detection of matrix size in e.g. matlab, python).

Implements oomph::Matrix< double, CRDoubleMatrix >.

     {
       CR_matrix.output_bottom_right_zero_helper(outfile);
     }

References CR_matrix, and oomph::CRMatrix< T >::output_bottom_right_zero_helper().

◆ redistribute()

void oomph::CRDoubleMatrix::redistribute ( const LinearAlgebraDistribution *const & dist_pt )

The contents of the matrix are redistributed to match the new distribution. In a non-MPI build this method does nothing. NOTE 1: The current distribution and the new distribution must have the same number of global rows. NOTE 2: The current distribution and the new distribution must have the same Communicator.

   {
 #ifdef OOMPH_HAS_MPI
 #ifdef PARANOID
     // paranoid check that the nrows for both distributions is the
     // same
     if (dist_pt->nrow() != this->distribution_pt()->nrow())
     {
       std::ostringstream error_message;
       error_message << "The number of global rows in the new distribution ("
                     << dist_pt->nrow() << ") is not equal to the number"
                     << " of global rows in the current distribution ("
                     << this->distribution_pt()->nrow() << ").\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
     // paranoid check that the current distribution and the new distribution
     // have the same Communicator
     OomphCommunicator temp_comm(*dist_pt->communicator_pt());
     if (!(temp_comm == *this->distribution_pt()->communicator_pt()))
     {
       std::ostringstream error_message;
       error_message << "The new distribution and the current distribution must "
                     << "have the same communicator.";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
     // paranoid check that the matrix is build
     if (!this->built())
     {
       std::ostringstream error_message;
       error_message << "The matrix must be build to be redistributed";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // if the two distributions are not the same
     // =========================================
     if (!((*this->distribution_pt()) == *dist_pt))
     {
       // Get the number of columns to build the matrix.
       unsigned long ncol = this->ncol();
  
       // current data
       int* current_row_start = this->row_start();
       int* current_column_index = this->column_index();
       double* current_value = this->value();
  
       // get the rank and the number of processors
       int my_rank = this->distribution_pt()->communicator_pt()->my_rank();
       int nproc = this->distribution_pt()->communicator_pt()->nproc();
  
       // if both distributions are distributed
       // =====================================
       if (this->distributed() && dist_pt->distributed())
       {
         // new nrow_local and first_row data
         Vector<unsigned> new_first_row(nproc);
         Vector<unsigned> new_nrow_local(nproc);
         Vector<unsigned> current_first_row(nproc);
         Vector<unsigned> current_nrow_local(nproc);
         for (int i = 0; i < nproc; i++)
         {
           new_first_row[i] = dist_pt->first_row(i);
           new_nrow_local[i] = dist_pt->nrow_local(i);
           current_first_row[i] = this->first_row(i);
           current_nrow_local[i] = this->nrow_local(i);
         }
  
         // compute which local rows are expected to be received from each
         // processor / sent to each processor
         Vector<unsigned> first_row_for_proc(nproc, 0);
         Vector<unsigned> nrow_local_for_proc(nproc, 0);
         Vector<unsigned> first_row_from_proc(nproc, 0);
         Vector<unsigned> nrow_local_from_proc(nproc, 0);
  
         // for every processor compute first_row and nrow_local that will
         // will sent and received by this processor
         for (int p = 0; p < nproc; p++)
         {
           // start with data to be sent
           if ((new_first_row[p] <
                (current_first_row[my_rank] + current_nrow_local[my_rank])) &&
               (current_first_row[my_rank] <
                (new_first_row[p] + new_nrow_local[p])))
           {
             first_row_for_proc[p] =
               std::max(current_first_row[my_rank], new_first_row[p]);
             nrow_local_for_proc[p] =
               std::min(
                 (current_first_row[my_rank] + current_nrow_local[my_rank]),
                 (new_first_row[p] + new_nrow_local[p])) -
               first_row_for_proc[p];
           }
  
           // and data to be received
           if ((new_first_row[my_rank] <
                (current_first_row[p] + current_nrow_local[p])) &&
               (current_first_row[p] <
                (new_first_row[my_rank] + new_nrow_local[my_rank])))
           {
             first_row_from_proc[p] =
               std::max(current_first_row[p], new_first_row[my_rank]);
             nrow_local_from_proc[p] =
               std::min((current_first_row[p] + current_nrow_local[p]),
                        (new_first_row[my_rank] + new_nrow_local[my_rank])) -
               first_row_from_proc[p];
           }
         }
  
         // determine how many nnzs to send to each processor
         Vector<unsigned> nnz_for_proc(nproc, 0);
         for (int p = 0; p < nproc; p++)
         {
           if (nrow_local_for_proc[p] > 0)
           {
             nnz_for_proc[p] = (current_row_start[first_row_for_proc[p] -
                                                  current_first_row[my_rank] +
                                                  nrow_local_for_proc[p]] -
                                current_row_start[first_row_for_proc[p] -
                                                  current_first_row[my_rank]]);
           }
         }
  
         // next post non-blocking sends and recv for the nnzs
         Vector<unsigned> nnz_from_proc(nproc, 0);
         Vector<MPI_Request> send_req;
         Vector<MPI_Request> nnz_recv_req;
         for (int p = 0; p < nproc; p++)
         {
           if (p != my_rank)
           {
             // send
             if (nrow_local_for_proc[p] > 0)
             {
               MPI_Request req;
               MPI_Isend(&nnz_for_proc[p],
                         1,
                         MPI_UNSIGNED,
                         p,
                         0,
                         this->distribution_pt()->communicator_pt()->mpi_comm(),
                         &req);
               send_req.push_back(req);
             }
  
             // recv
             if (nrow_local_from_proc[p] > 0)
             {
               MPI_Request req;
               MPI_Irecv(&nnz_from_proc[p],
                         1,
                         MPI_UNSIGNED,
                         p,
                         0,
                         this->distribution_pt()->communicator_pt()->mpi_comm(),
                         &req);
               nnz_recv_req.push_back(req);
             }
           }
           // "send to self"
           else
           {
             nnz_from_proc[p] = nnz_for_proc[p];
           }
         }
  
         // allocate new storage for the new row_start
         int* new_row_start = new int[new_nrow_local[my_rank] + 1];
  
         // wait for recvs to complete
         unsigned n_recv_req = nnz_recv_req.size();
         if (n_recv_req > 0)
         {
           Vector<MPI_Status> recv_status(n_recv_req);
           MPI_Waitall(n_recv_req, &nnz_recv_req[0], &recv_status[0]);
         }
  
         // compute the nnz offset for each processor
         unsigned next_row = 0;
         unsigned nnz_count = 0;
         Vector<unsigned> nnz_offset(nproc, 0);
         for (int p = 0; p < nproc; p++)
         {
           unsigned pp = 0;
           while (new_first_row[pp] != next_row)
           {
             pp++;
           }
           nnz_offset[pp] = nnz_count;
           nnz_count += nnz_from_proc[pp];
           next_row += new_nrow_local[pp];
         }
  
         // allocate storage for the values and column indices
         int* new_column_index = new int[nnz_count];
         double* new_value = new double[nnz_count];
  
         // post the sends and recvs for the matrix data
         Vector<MPI_Request> recv_req;
         MPI_Aint base_address;
         MPI_Get_address(new_value, &base_address);
         for (int p = 0; p < nproc; p++)
         {
           // communicated with other processors
           if (p != my_rank)
           {
             // SEND
             if (nrow_local_for_proc[p] > 0)
             {
               // array of datatypes
               MPI_Datatype types[3];
  
               // array of offsets
               MPI_Aint offsets[3];
  
               // array of lengths
               int len[3];
  
               // row start
               unsigned first_row_to_send =
                 first_row_for_proc[p] - current_first_row[my_rank];
               MPI_Type_contiguous(nrow_local_for_proc[p], MPI_INT, &types[0]);
               MPI_Type_commit(&types[0]);
               len[0] = 1;
               MPI_Get_address(current_row_start + first_row_to_send,
                               &offsets[0]);
               offsets[0] -= base_address;
  
               // values
               unsigned first_coef_to_send =
                 current_row_start[first_row_to_send];
               MPI_Type_contiguous(nnz_for_proc[p], MPI_DOUBLE, &types[1]);
               MPI_Type_commit(&types[1]);
               len[1] = 1;
               MPI_Get_address(current_value + first_coef_to_send, &offsets[1]);
               offsets[1] -= base_address;
  
               // column index
               MPI_Type_contiguous(nnz_for_proc[p], MPI_DOUBLE, &types[2]);
               MPI_Type_commit(&types[2]);
               len[2] = 1;
               MPI_Get_address(current_column_index + first_coef_to_send,
                               &offsets[2]);
               offsets[2] -= base_address;
  
               // build the combined datatype
               MPI_Datatype send_type;
               MPI_Type_create_struct(3, len, offsets, types, &send_type);
               MPI_Type_commit(&send_type);
               MPI_Type_free(&types[0]);
               MPI_Type_free(&types[1]);
               MPI_Type_free(&types[2]);
  
               // and send
               MPI_Request req;
               MPI_Isend(new_value,
                         1,
                         send_type,
                         p,
                         1,
                         this->distribution_pt()->communicator_pt()->mpi_comm(),
                         &req);
               send_req.push_back(req);
               MPI_Type_free(&send_type);
             }
  
             // RECV
             if (nrow_local_from_proc[p] > 0)
             {
               // array of datatypes
               MPI_Datatype types[3];
  
               // array of offsets
               MPI_Aint offsets[3];
  
               // array of lengths
               int len[3];
  
               // row start
               unsigned first_row_to_recv =
                 first_row_from_proc[p] - new_first_row[my_rank];
               MPI_Type_contiguous(nrow_local_from_proc[p], MPI_INT, &types[0]);
               MPI_Type_commit(&types[0]);
               len[0] = 1;
               MPI_Get_address(new_row_start + first_row_to_recv, &offsets[0]);
               offsets[0] -= base_address;
  
               // values
               unsigned first_coef_to_recv = nnz_offset[p];
               MPI_Type_contiguous(nnz_from_proc[p], MPI_DOUBLE, &types[1]);
               MPI_Type_commit(&types[1]);
               len[1] = 1;
               MPI_Get_address(new_value + first_coef_to_recv, &offsets[1]);
               offsets[1] -= base_address;
  
               // column index
               MPI_Type_contiguous(nnz_from_proc[p], MPI_INT, &types[2]);
               MPI_Type_commit(&types[2]);
               len[2] = 1;
               MPI_Get_address(new_column_index + first_coef_to_recv,
                               &offsets[2]);
               offsets[2] -= base_address;
  
               // build the combined datatype
               MPI_Datatype recv_type;
               MPI_Type_create_struct(3, len, offsets, types, &recv_type);
               MPI_Type_commit(&recv_type);
               MPI_Type_free(&types[0]);
               MPI_Type_free(&types[1]);
               MPI_Type_free(&types[2]);
  
               // and send
               MPI_Request req;
               MPI_Irecv(new_value,
                         1,
                         recv_type,
                         p,
                         1,
                         this->distribution_pt()->communicator_pt()->mpi_comm(),
                         &req);
               recv_req.push_back(req);
               MPI_Type_free(&recv_type);
             }
           }
           // other wise transfer data internally
           else
           {
             unsigned j =
               first_row_for_proc[my_rank] - current_first_row[my_rank];
             unsigned k = first_row_from_proc[my_rank] - new_first_row[my_rank];
             for (unsigned i = 0; i < nrow_local_for_proc[my_rank]; i++)
             {
               new_row_start[k + i] = current_row_start[j + i];
             }
             unsigned first_coef_to_send = current_row_start[j];
             for (unsigned i = 0; i < nnz_for_proc[my_rank]; i++)
             {
               new_value[nnz_offset[p] + i] =
                 current_value[first_coef_to_send + i];
               new_column_index[nnz_offset[p] + i] =
                 current_column_index[first_coef_to_send + i];
             }
           }
         }
  
         // wait for all recvs to complete
         n_recv_req = recv_req.size();
         if (n_recv_req > 0)
         {
           Vector<MPI_Status> recv_status(n_recv_req);
           MPI_Waitall(n_recv_req, &recv_req[0], &recv_status[0]);
         }
  
         // next we need to update the row starts
         for (int p = 0; p < nproc; p++)
         {
           if (nrow_local_from_proc[p] > 0)
           {
             unsigned first_row =
               first_row_from_proc[p] - new_first_row[my_rank];
             unsigned last_row = first_row + nrow_local_from_proc[p] - 1;
             int update = nnz_offset[p] - new_row_start[first_row];
             for (unsigned i = first_row; i <= last_row; i++)
             {
               new_row_start[i] += update;
             }
           }
         }
         new_row_start[dist_pt->nrow_local()] = nnz_count;
  
         // wait for sends to complete
         unsigned n_send_req = send_req.size();
         if (n_recv_req > 0)
         {
           Vector<MPI_Status> send_status(n_recv_req);
           MPI_Waitall(n_send_req, &send_req[0], &send_status[0]);
         }
         // if (my_rank == 0)
         //  {
         //   CRDoubleMatrix* m_pt = this->global_matrix();
         //   m_pt->sparse_indexed_output("m1.dat");
         //  }
  
         //
         this->build(dist_pt);
         this->build_without_copy(
           ncol, nnz_count, new_value, new_column_index, new_row_start);
         // if (my_rank == 0)
         //  {
         //   CRDoubleMatrix* m_pt = this->global_matrix();
         //   m_pt->sparse_indexed_output("m2.dat");
         //  }
         //       this->sparse_indexed_output(oomph_info);
         abort();
       }
  
       // if this matrix is distributed but the new distributed matrix is global
       // ======================================================================
       else if (this->distributed() && !dist_pt->distributed())
       {
         // nnz
         int nnz = this->nnz();
  
         // nrow global
         unsigned nrow = this->nrow();
  
         // cache nproc
         int nproc = this->distribution_pt()->communicator_pt()->nproc();
  
         // get the nnzs on the other processors
         int* dist_nnz_pt = new int[nproc];
         MPI_Allgather(&nnz,
                       1,
                       MPI_INT,
                       dist_nnz_pt,
                       1,
                       MPI_INT,
                       this->distribution_pt()->communicator_pt()->mpi_comm());
  
         // create an int array of first rows and nrow local and
         // compute nnz global
         int* dist_first_row = new int[nproc];
         int* dist_nrow_local = new int[nproc];
         for (int p = 0; p < nproc; p++)
         {
           dist_first_row[p] = this->first_row(p);
           dist_nrow_local[p] = this->nrow_local(p);
         }
  
         // conpute the offset for the values and column index data
         // compute the nnz offset for each processor
         int next_row = 0;
         unsigned nnz_count = 0;
         Vector<unsigned> nnz_offset(nproc, 0);
         for (int p = 0; p < nproc; p++)
         {
           unsigned pp = 0;
           while (dist_first_row[pp] != next_row)
           {
             pp++;
           }
           nnz_offset[pp] = nnz_count;
           nnz_count += dist_nnz_pt[pp];
           next_row += dist_nrow_local[pp];
         }
  
         // get pointers to the (current) distributed data
         int* dist_row_start = this->row_start();
         int* dist_column_index = this->column_index();
         double* dist_value = this->value();
  
         // space for the global matrix
         int* global_row_start = new int[nrow + 1];
         int* global_column_index = new int[nnz_count];
         double* global_value = new double[nnz_count];
  
         // post the sends and recvs for the matrix data
         Vector<MPI_Request> recv_req;
         Vector<MPI_Request> send_req;
         MPI_Aint base_address;
         MPI_Get_address(global_value, &base_address);
  
         // SEND
         if (dist_nrow_local[my_rank] > 0)
         {
           // types
           MPI_Datatype types[3];
  
           // offsets
           MPI_Aint offsets[3];
  
           // lengths
           int len[3];
  
           // row start
           MPI_Type_contiguous(dist_nrow_local[my_rank], MPI_INT, &types[0]);
           MPI_Type_commit(&types[0]);
           MPI_Get_address(dist_row_start, &offsets[0]);
           offsets[0] -= base_address;
           len[0] = 1;
  
           // value
           MPI_Type_contiguous(nnz, MPI_DOUBLE, &types[1]);
           MPI_Type_commit(&types[1]);
           MPI_Get_address(dist_value, &offsets[1]);
           offsets[1] -= base_address;
           len[1] = 1;
  
           // column indices
           MPI_Type_contiguous(nnz, MPI_INT, &types[2]);
           MPI_Type_commit(&types[2]);
           MPI_Get_address(dist_column_index, &offsets[2]);
           offsets[2] -= base_address;
           len[2] = 1;
  
           // build the send type
           MPI_Datatype send_type;
           MPI_Type_create_struct(3, len, offsets, types, &send_type);
           MPI_Type_commit(&send_type);
           MPI_Type_free(&types[0]);
           MPI_Type_free(&types[1]);
           MPI_Type_free(&types[2]);
  
           // and send
           for (int p = 0; p < nproc; p++)
           {
             if (p != my_rank)
             {
               MPI_Request req;
               MPI_Isend(global_value,
                         1,
                         send_type,
                         p,
                         1,
                         this->distribution_pt()->communicator_pt()->mpi_comm(),
                         &req);
               send_req.push_back(req);
             }
           }
           MPI_Type_free(&send_type);
         }
  
         // RECV
         for (int p = 0; p < nproc; p++)
         {
           // communicated with other processors
           if (p != my_rank)
           {
             // RECV
             if (dist_nrow_local[p] > 0)
             {
               // types
               MPI_Datatype types[3];
  
               // offsets
               MPI_Aint offsets[3];
  
               // lengths
               int len[3];
  
               // row start
               MPI_Type_contiguous(dist_nrow_local[p], MPI_INT, &types[0]);
               MPI_Type_commit(&types[0]);
               MPI_Get_address(global_row_start + dist_first_row[p],
                               &offsets[0]);
               offsets[0] -= base_address;
               len[0] = 1;
  
               // value
               MPI_Type_contiguous(dist_nnz_pt[p], MPI_DOUBLE, &types[1]);
               MPI_Type_commit(&types[1]);
               MPI_Get_address(global_value + nnz_offset[p], &offsets[1]);
               offsets[1] -= base_address;
               len[1] = 1;
  
               // column indices
               MPI_Type_contiguous(dist_nnz_pt[p], MPI_INT, &types[2]);
               MPI_Type_commit(&types[2]);
               MPI_Get_address(global_column_index + nnz_offset[p], &offsets[2]);
               offsets[2] -= base_address;
               len[2] = 1;
  
               // build the send type
               MPI_Datatype recv_type;
               MPI_Type_create_struct(3, len, offsets, types, &recv_type);
               MPI_Type_commit(&recv_type);
               MPI_Type_free(&types[0]);
               MPI_Type_free(&types[1]);
               MPI_Type_free(&types[2]);
  
               // and send
               MPI_Request req;
               MPI_Irecv(global_value,
                         1,
                         recv_type,
                         p,
                         1,
                         this->distribution_pt()->communicator_pt()->mpi_comm(),
                         &req);
               recv_req.push_back(req);
               MPI_Type_free(&recv_type);
             }
           }
           // otherwise send to self
           else
           {
             unsigned nrow_local = dist_nrow_local[my_rank];
             unsigned first_row = dist_first_row[my_rank];
             for (unsigned i = 0; i < nrow_local; i++)
             {
               global_row_start[first_row + i] = dist_row_start[i];
             }
             unsigned offset = nnz_offset[my_rank];
             for (int i = 0; i < nnz; i++)
             {
               global_value[offset + i] = dist_value[i];
               global_column_index[offset + i] = dist_column_index[i];
             }
           }
         }
  
         // wait for all recvs to complete
         unsigned n_recv_req = recv_req.size();
         if (n_recv_req > 0)
         {
           Vector<MPI_Status> recv_status(n_recv_req);
           MPI_Waitall(n_recv_req, &recv_req[0], &recv_status[0]);
         }
  
         // finally the last row start
         global_row_start[nrow] = nnz_count;
  
         // update the other row start
         for (int p = 0; p < nproc; p++)
         {
           for (int i = 0; i < dist_nrow_local[p]; i++)
           {
             unsigned j = dist_first_row[p] + i;
             global_row_start[j] += nnz_offset[p];
           }
         }
  
         // wait for sends to complete
         unsigned n_send_req = send_req.size();
         if (n_recv_req > 0)
         {
           Vector<MPI_Status> send_status(n_recv_req);
           MPI_Waitall(n_send_req, &send_req[0], &send_status[0]);
         }
  
         // rebuild the matrix
         LinearAlgebraDistribution* dist_pt = new LinearAlgebraDistribution(
           this->distribution_pt()->communicator_pt(), nrow, false);
         this->build(dist_pt);
         this->build_without_copy(
           ncol, nnz_count, global_value, global_column_index, global_row_start);
  
         // clean up
         delete dist_pt;
         delete[] dist_first_row;
         delete[] dist_nrow_local;
         delete[] dist_nnz_pt;
       }
  
       // other the matrix is not distributed but it needs to be turned
       // into a distributed matrix
       // =============================================================
       else if (!this->distributed() && dist_pt->distributed())
       {
         // cache the new nrow_local
         unsigned nrow_local = dist_pt->nrow_local();
  
         // and first_row
         unsigned first_row = dist_pt->first_row();
  
         // get pointers to the (current) distributed data
         int* global_row_start = this->row_start();
         int* global_column_index = this->column_index();
         double* global_value = this->value();
  
         // determine the number of non zeros required by this processor
         unsigned nnz = global_row_start[first_row + nrow_local] -
                        global_row_start[first_row];
  
         // allocate
         int* dist_row_start = new int[nrow_local + 1];
         int* dist_column_index = new int[nnz];
         double* dist_value = new double[nnz];
  
         // copy
         int offset = global_row_start[first_row];
         for (unsigned i = 0; i <= nrow_local; i++)
         {
           dist_row_start[i] = global_row_start[first_row + i] - offset;
         }
         for (unsigned i = 0; i < nnz; i++)
         {
           dist_column_index[i] = global_column_index[offset + i];
           dist_value[i] = global_value[offset + i];
         }
  
         // rebuild
         this->build(dist_pt);
         this->build_without_copy(
           ncol, nnz, dist_value, dist_column_index, dist_row_start);
       }
     }
 #endif
   }

References build(), build_without_copy(), built(), column_index(), oomph::LinearAlgebraDistribution::communicator_pt(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearAlgebraDistribution::first_row(), oomph::DistributableLinearAlgebraObject::first_row(), i, j, k, max, min, oomph::OomphCommunicator::my_rank(), ncol(), nnz(), oomph::OomphCommunicator::nproc(), oomph::LinearAlgebraDistribution::nrow(), nrow(), oomph::LinearAlgebraDistribution::nrow_local(), oomph::DistributableLinearAlgebraObject::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p, row_start(), and value().

Referenced by oomph::CRDoubleMatrixHelpers::create_uniformly_distributed_matrix(), oomph::Problem::get_eigenproblem_matrices(), and oomph::Problem::get_jacobian().

◆ row_start() [1/2]

int* oomph::CRDoubleMatrix::row_start ( )

inline

Access to C-style row_start array.

     {
       return CR_matrix.row_start();
     }

References CR_matrix, and oomph::CRMatrix< T >::row_start().

Referenced by add(), build(), build_without_copy(), oomph::CRDoubleMatrixHelpers::concatenate(), CRDoubleMatrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix_for_aztecoo(), oomph::HypreHelpers::create_HYPRE_Matrix(), oomph::CRDoubleMatrixHelpers::deep_copy(), entries_are_sorted(), oomph::MumpsSolver::factorise(), oomph::SuperLUSolver::factorise_serial(), oomph::CRDoubleMatrixHelpers::gershgorin_eigenvalue_estimate(), get_matrix_transpose(), global_matrix(), inf_norm(), oomph::CRDoubleMatrixHelpers::inf_norm(), oomph::BlockPreconditioner< MATRIX >::internal_get_block(), matrix_reduction(), oomph::MGSolver< DIM >::modify_restriction_matrices(), multiply(), multiply_transpose(), redistribute(), oomph::MatrixBasedLumpedPreconditioner< MATRIX >::setup(), oomph::ILUZeroPreconditioner< CRDoubleMatrix >::setup(), oomph::PseudoElasticPreconditionerSubsidiaryPreconditionerOld::setup(), oomph::HelmholtzMGPreconditioner< DIM >::setup_coarsest_level_structures(), oomph::GS< CRDoubleMatrix >::solve_helper(), sort_entries(), and sparse_indexed_output_with_offset().

◆ row_start() [2/2]

const int* oomph::CRDoubleMatrix::row_start ( ) const

inline

Access to C-style row_start array (const version)

     {
       return CR_matrix.row_start();
     }

References CR_matrix, and oomph::CRMatrix< T >::row_start().

◆ serial_matrix_matrix_multiply_method() [1/2]

unsigned& oomph::CRDoubleMatrix::serial_matrix_matrix_multiply_method ( )

inline

Access function to Serial_matrix_matrix_multiply_method, the flag which determines the matrix matrix multiplication method used for serial matrices. Method 1: First runs through this matrix and matrix_in to find the storage requirements for result - arrays of the correct size are then allocated before performing the calculation. Minimises memory requirements but more costly. Method 2: Grows storage for values and column indices of result 'on the fly' using an array of maps. Faster but more memory intensive. Method 3: Grows storage for values and column indices of result 'on the fly' using a vector of vectors. Not particularly impressive on the platforms we tried... Method 4: Trilinos Epetra Matrix Matrix multiply. Method 5: Trilinos Epetra Matrix Matrix multiply (ML based).

     {
       return Serial_matrix_matrix_multiply_method;
     }

References Serial_matrix_matrix_multiply_method.

Referenced by oomph::CRDoubleMatrixHelpers::deep_copy(), and main().

◆ serial_matrix_matrix_multiply_method() [2/2]

const unsigned& oomph::CRDoubleMatrix::serial_matrix_matrix_multiply_method ( ) const

inline

Read only access function (const version) to Serial_matrix_matrix_multiply_method, the flag which determines the matrix matrix multiplication method used for serial matrices. Method 1: First runs through this matrix and matrix_in to find the storage requirements for result - arrays of the correct size are then allocated before performing the calculation. Minimises memory requirements but more costly. Method 2: Grows storage for values and column indices of result 'on the fly' using an array of maps. Faster but more memory intensive. Method 3: Grows storage for values and column indices of result 'on the fly' using a vector of vectors. Not particularly impressive on the platforms we tried... Method 4: Trilinos Epetra Matrix Matrix multiply. Method 5: Trilinos Epetra Matrix Matrix multiply (ML based).

     {
       return Serial_matrix_matrix_multiply_method;
     }

References Serial_matrix_matrix_multiply_method.

◆ sort_entries()

void oomph::CRDoubleMatrix::sort_entries ( )

Sorts the entries associated with each row of the matrix in the column index vector and the value vector into ascending order and sets up the Index_of_diagonal_entries vector

This helper function sorts the entries in the column index vector and the value vector. During the construction of the matrix the entries were most likely assigned in an arbitrary order. As a result, it cannot be assumed that the entries in the column index vector corresponding to each row of the matrix have been arranged in increasing order. During the setup an additional vector will be set up; Index_of_diagonal_entries. The i-th entry of this vector contains the index of the last entry below or on the diagonal. If there are no entries below or on the diagonal then the corresponding entry is -1. If, however, there are no entries in the row then the entry is irrelevant and is kept as the initialised value; 0.

   {
 #ifdef OOMPH_HAS_MPI
     // We only need to produce a warning if the matrix is distributed
     if (this->distributed())
     {
       // Create an ostringstream object to store the warning message
       std::ostringstream warning_message;
  
       // Create the warning messsage
       warning_message << "This method currently works for serial but "
                       << "has not been tested with MPI!\n";
  
       // Issue the warning
       OomphLibWarning(warning_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Get the number of rows in the matrix
     unsigned n_rows = this->nrow();
  
     // Acquire access to the value, row_start and column_index arrays from
     // the CR matrix. Since we do not change anything in row_start_pt we
     // give it the const prefix
     double* value_pt = this->value();
     int* column_index_pt = this->column_index();
     const int* row_start_pt = this->row_start();
  
     // Resize the Index_of_diagonal_entries vector
     Index_of_diagonal_entries.resize(n_rows, 0);
  
     // Vector of pairs to store the column_index of each value in the i-th row
     // and its corresponding matrix entry
     Vector<std::pair<int, double>> column_index_and_value_row_i;
  
     // Loop over the rows of the matrix
     for (unsigned i = 0; i < n_rows; i++)
     {
       // Find the number of nonzeros in the i-th row
       unsigned nnz_row_i = *(row_start_pt + i + 1) - *(row_start_pt + i);
  
       // Variable to store the start of the i-th row
       unsigned i_row_start = *(row_start_pt + i);
  
       // If there are no nonzeros in this row then the i-th entry of the vector
       // Index_of_diagonal_entries is irrelevant so we can simply let it be 0
       if (nnz_row_i == 0)
       {
         // Set the i-th entry
         Index_of_diagonal_entries[i] = 0;
       }
       // If there are nonzeros in the i-th row
       else
       {
         // If there is more than one entry in the row resize the vector
         // column_index_and_value_row_i
         column_index_and_value_row_i.resize(nnz_row_i);
  
         // Loop over the entries in the row
         for (unsigned j = 0; j < nnz_row_i; j++)
         {
           // Assign the appropriate entries to column_index_and_value_row_i
           column_index_and_value_row_i[j] =
             std::make_pair(*(column_index_pt + i_row_start + j),
                            *(value_pt + i_row_start + j));
         }
  
         // Sort the vector of pairs using the struct
         // CRDoubleMatrixComparisonHelper
         std::sort(column_index_and_value_row_i.begin(),
                   column_index_and_value_row_i.end(),
                   Comparison_struct);
  
         //-----------------------------------------------------------------------
         // Now that the entries of the i-th row have been sorted we can read
         // them back into value_pt and column_index_pt:
         //-----------------------------------------------------------------------
  
         // Create a boolean variable to indicate whether or not the i-th entry
         // of Index_of_diagonal_entries has been set
         bool is_ith_entry_set = false;
  
         // Loop over the entries in the vector column_index_and_value_row_i and
         // assign its entries to value_pt and column_index_pt
         for (unsigned j = 0; j < nnz_row_i; j++)
         {
           // Set the column index of the j-th nonzero value in the i-th row of
           // the matrix
           *(column_index_pt + i_row_start + j) =
             column_index_and_value_row_i[j].first;
  
           // Set the value of the j-th nonzero value in the i-th row of
           // the matrix
           *(value_pt + i_row_start + j) =
             column_index_and_value_row_i[j].second;
  
           // This if statement is used to set the i-th entry of the vector
           // Index_of_diagonal_entries if it has not yet been set
           if (!is_ith_entry_set)
           {
             // If the column index of the first entry in row i is greater than
             // the row number then the first entry must lie above the diagonal
             if (unsigned(*(column_index_pt + i_row_start)) > i)
             {
               // If the column index of the first entry in the row is greater
               // than the row number, i, then the i-th entry of
               // Index_of_diagonal_entries needs to be set to -1 to indicate
               // there are no entries below or on the diagonal
               Index_of_diagonal_entries[i] = -1;
  
               // Indicate that the i-th entry of Index_of_diagonal_entries has
               // been set
               is_ith_entry_set = true;
             }
             // If there are entries below or on the diagonal
             else
             {
               // If there is only one entry in the row then we know that this
               // will be the last entry below or on the diagonal because we have
               // eliminated the possibility that if there is only one entry,
               // that it lies above the diagonal
               if (nnz_row_i == 1)
               {
                 // Set the index of the current entry to be the value of i-th
                 // entry of Index_of_diagonal_entries
                 Index_of_diagonal_entries[i] = i_row_start + j;
  
                 // Indicate that the i-th entry of Index_of_diagonal_entries has
                 // been set
                 is_ith_entry_set = true;
               }
               // It remains to now check the case that there is more than one
               // entry in the row. If there is more than one entry in the row
               // and there are entries below or on the diagonal then we have
               // three cases:
               //          (1) The current entry lies on the diagonal;
               //          (2) The current entry lies above the diagonal;
               //          (3) The current entry lies below the diagonal;
               // The first case can easily be checked as done below. If the
               // second case occurs then we have just passed the last entry. We
               // know this because at least one entry lies on or below the
               // diagonal. If the second case it true then we need to assign the
               // previous entry to the vector Index_of_diagonal_entries.
               // Finally, we are left with case (3), which can be split into two
               // cases:
               //          (3.1) The current entry lies below the diagonal but it
               //                is not the last entry below or on the diagonal;
               //          (3.2) The current entry lies below the diagonal and is
               //                the last entry below or on the diagonal.
               // If case (3.1) holds then we can simply wait until we get to the
               // next entry in the row and examine that. If the next entry lies
               // on the diagonal then it will be swept up by case (1). If the
               // next entry lies above the diagonal then case (2) will sweep it
               // up and if neither is the case then we wait until the next entry
               // and so on. If, instead, case (3.2) holds then our last check
               // simply needs to check if the current entry is the last entry in
               // the row because if the last entry lies on the diagonal, case
               // (1) will sweep it up. If it lies above the diagonal, case (2)
               // will take care of it. Therefore, the only remaining case is
               // that it lies strictly below the diagonal and since it is the
               // last entry in the row it means the index of this entry needs to
               // be assigned to Index_of_diagonal_entries
  
               // Case (1) : The current entry lies on the diagonal
               else if (unsigned(*(column_index_pt + i_row_start + j)) == i)
               {
                 // Set the index of the current entry to be the value of i-th
                 // entry of Index_of_diagonal_entries
                 Index_of_diagonal_entries[i] = i_row_start + j;
  
                 // Indicate that the i-th entry of Index_of_diagonal_entries has
                 // been set
                 is_ith_entry_set = true;
               }
               // Case (2) : The current entry lies above the diagonal
               else if (unsigned(*(column_index_pt + i_row_start + j)) > i)
               {
                 // Set the index of the current entry to be the value of i-th
                 // entry of Index_of_diagonal_entries
                 Index_of_diagonal_entries[i] = i_row_start + j - 1;
  
                 // Indicate that the i-th entry of Index_of_diagonal_entries has
                 // been set
                 is_ith_entry_set = true;
               }
               // Case (3.2) : The current entry is the last entry in the row
               else if (j == nnz_row_i - 1)
               {
                 // Set the index of the current entry to be the value of i-th
                 // entry of Index_of_diagonal_entries
                 Index_of_diagonal_entries[i] = i_row_start + j;
  
                 // Indicate that the i-th entry of Index_of_diagonal_entries has
                 // been set
                 is_ith_entry_set = true;
               } // if (nnz_row_i==1) else if
             } // if (*(column_index_pt+i_row_start)>i)
           } // if (!is_ith_entry_set)
         } // for (unsigned j=0;j<nnz_row_i;j++)
       } // if (nnz_row_i==0) else
     } // for (unsigned i=0;i<n_rows;i++)
   } // End of sort_entries()

References column_index(), Comparison_struct, oomph::DistributableLinearAlgebraObject::distributed(), i, Index_of_diagonal_entries, j, nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, row_start(), and value().

◆ sparse_indexed_output_helper()

void oomph::CRDoubleMatrix::sparse_indexed_output_helper ( std::ostream & outfile ) const

inlinevirtual

Indexed output function to print a matrix to the stream outfile as i,j,a(i,j) for a(i,j)!=0 only.

Implements oomph::Matrix< double, CRDoubleMatrix >.

     {
       CR_matrix.sparse_indexed_output_helper(outfile);
     }

References CR_matrix, and oomph::CRMatrix< T >::sparse_indexed_output_helper().

◆ sparse_indexed_output_with_offset()

void oomph::CRDoubleMatrix::sparse_indexed_output_with_offset ( std::string filename )

inline

Indexed output function to print a matrix to a file as i,j,a(i,j) for a(i,j)!=0 only. Specify filename. This uses acual global row numbers.

     {
       // Get offset
       unsigned first_row = distribution_pt()->first_row();
  
       // Open file
       std::ofstream some_file;
       some_file.open(filename.c_str());
       unsigned n = nrow_local();
       for (unsigned long i = 0; i < n; i++)
       {
         for (long j = row_start()[i]; j < row_start()[i + 1]; j++)
         {
           some_file << first_row + i << " " << column_index()[j] << " "
                     << value()[j] << std::endl;
         }
       }
       some_file.close();
     }

References column_index(), oomph::DistributableLinearAlgebraObject::distribution_pt(), MergeRestartFiles::filename, oomph::LinearAlgebraDistribution::first_row(), oomph::DistributableLinearAlgebraObject::first_row(), i, j, n, oomph::DistributableLinearAlgebraObject::nrow_local(), row_start(), and value().

Referenced by oomph::SpaceTimeNavierStokesSubsidiaryPreconditioner::setup().

◆ value() [1/2]

double* oomph::CRDoubleMatrix::value ( )

inline

Access to C-style value array.

     {
       return CR_matrix.value();
     }

References CR_matrix, and oomph::SparseMatrix< T, MATRIX_TYPE >::value().

Referenced by add(), build(), build_without_copy(), oomph::CRDoubleMatrixHelpers::concatenate(), CRDoubleMatrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix(), oomph::TrilinosEpetraHelpers::create_distributed_epetra_matrix_for_aztecoo(), oomph::HypreHelpers::create_HYPRE_Matrix(), oomph::CRDoubleMatrixHelpers::deep_copy(), oomph::MumpsSolver::factorise(), oomph::SuperLUSolver::factorise_serial(), oomph::CRDoubleMatrixHelpers::gershgorin_eigenvalue_estimate(), get_matrix_transpose(), global_matrix(), inf_norm(), oomph::CRDoubleMatrixHelpers::inf_norm(), oomph::BlockPreconditioner< MATRIX >::internal_get_block(), matrix_reduction(), oomph::MGSolver< DIM >::modify_restriction_matrices(), multiply(), multiply_transpose(), redistribute(), oomph::MatrixBasedLumpedPreconditioner< MATRIX >::setup(), oomph::ILUZeroPreconditioner< CRDoubleMatrix >::setup(), oomph::PseudoElasticPreconditionerSubsidiaryPreconditionerOld::setup(), oomph::HelmholtzMGPreconditioner< DIM >::setup_coarsest_level_structures(), oomph::GS< CRDoubleMatrix >::solve_helper(), sort_entries(), and sparse_indexed_output_with_offset().

◆ value() [2/2]

const double* oomph::CRDoubleMatrix::value ( ) const

inline

Access to C-style value array (const version)

     {
       return CR_matrix.value();
     }

References CR_matrix, and oomph::SparseMatrix< T, MATRIX_TYPE >::value().

Member Data Documentation

◆ Built

bool oomph::CRDoubleMatrix::Built

private

Flag to indicate whether the matrix has been built - i.e. the distribution has been setup AND the matrix has been assembled.

Referenced by build(), build_without_copy(), built(), clear(), CRDoubleMatrix(), ludecompose(), multiply(), and multiply_transpose().

◆ Comparison_struct

struct oomph::CRDoubleMatrix::CRDoubleMatrixComparisonHelper oomph::CRDoubleMatrix::Comparison_struct

Referenced by sort_entries().

◆ CR_matrix

CRMatrix<double> oomph::CRDoubleMatrix::CR_matrix

private

Storage for the Matrix in CR Format.

Referenced by build(), build_without_copy(), clear(), column_index(), CRDoubleMatrix(), diagonal_entries(), inf_norm(), matrix_reduction(), multiply(), multiply_transpose(), ncol(), nnz(), operator()(), output_bottom_right_zero_helper(), row_start(), sparse_indexed_output_helper(), and value().

◆ Distributed_matrix_matrix_multiply_method

unsigned oomph::CRDoubleMatrix::Distributed_matrix_matrix_multiply_method

private

Flag to determine which matrix-matrix multiplication method is used (for distributed matrices)

Referenced by distributed_matrix_matrix_multiply_method().

◆ Index_of_diagonal_entries

Vector<int> oomph::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

Referenced by get_index_of_diagonal_entries(), and sort_entries().

◆ Serial_matrix_matrix_multiply_method

unsigned oomph::CRDoubleMatrix::Serial_matrix_matrix_multiply_method

private

Flag to determine which matrix-matrix multiplication method is used (for serial (or global) matrices)

Referenced by CRDoubleMatrix(), multiply(), and serial_matrix_matrix_multiply_method().

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

Classes

Public Member Functions

Public Attributes

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ CRDoubleMatrix() [1/4]

◆ CRDoubleMatrix() [2/4]

◆ CRDoubleMatrix() [3/4]

◆ CRDoubleMatrix() [4/4]

◆ ~CRDoubleMatrix()

Member Function Documentation

◆ add()

◆ build() [1/3]

◆ build() [2/3]

◆ build() [3/3]

◆ build_without_copy()

◆ built()

◆ clear()

◆ column_index() [1/2]

◆ column_index() [2/2]

◆ diagonal_entries()

◆ distributed_matrix_matrix_multiply_method() [1/2]

◆ distributed_matrix_matrix_multiply_method() [2/2]

◆ entries_are_sorted()

◆ get_index_of_diagonal_entries()

◆ get_matrix_transpose()

◆ global_matrix()

◆ inf_norm()

◆ lubksub()

◆ ludecompose()

◆ matrix_reduction()

◆ multiply() [1/2]

◆ multiply() [2/2]

◆ multiply_transpose()

◆ ncol()

◆ nnz()

◆ nrow()

◆ operator()()

◆ operator=()

◆ output_bottom_right_zero_helper()

◆ redistribute()

◆ row_start() [1/2]

◆ row_start() [2/2]

◆ serial_matrix_matrix_multiply_method() [1/2]

◆ serial_matrix_matrix_multiply_method() [2/2]

◆ sort_entries()

◆ sparse_indexed_output_helper()

◆ sparse_indexed_output_with_offset()

◆ value() [1/2]

◆ value() [2/2]

Member Data Documentation

◆ Built

◆ Comparison_struct

◆ CR_matrix

◆ Distributed_matrix_matrix_multiply_method

◆ Index_of_diagonal_entries

◆ Serial_matrix_matrix_multiply_method