A class for compressed row matrices. More...

#include <complex_matrices.h>

Inheritance diagram for oomph::CRComplexMatrix:

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

	CRComplexMatrix (const Vector< std::complex< double >> &value, const Vector< int > &column_index, const Vector< int > &row_start, const unsigned long &n, const unsigned long &m)

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

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

virtual	~CRComplexMatrix ()
	Destructor: Kill the LU decomposition if it has been computed. More...

void	enable_doc_stats ()

void	disable_doc_stats ()
	the solve 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...

std::complex< double >	operator() (const unsigned long &i, const unsigned long &j) const
	Overload the round-bracket access operator for read-only access. More...

int	ludecompose ()
	LU decomposition using SuperLU. More...

void	lubksub (Vector< std::complex< double >> &rhs)
	LU back solve for given RHS. More...

void	clean_up_memory ()
	LU clean up (perhaps this should happen in the destructor) More...

void	residual (const Vector< std::complex< double >> &x, const Vector< std::complex< double >> &b, Vector< std::complex< double >> &residual)
	Find the residual to x of Ax=b, i.e. r=b-Ax. More...

void	multiply (const Vector< std::complex< double >> &x, Vector< std::complex< double >> &soln)
	Multiply the matrix by the vector x: soln=Ax. More...

void	multiply_transpose (const Vector< std::complex< double >> &x, Vector< std::complex< double >> &soln)
	Multiply the transposed matrix by the vector x: soln=A^T x. More...

Public Member Functions inherited from oomph::CRMatrix< std::complex< double > >
	CRMatrix ()
	Default constructor. More...

	CRMatrix (const Vector< std::complex< double > > &value, const Vector< int > &column_index_, const Vector< int > &row_start_, const unsigned long &n, const unsigned long &m)

	CRMatrix (const CRMatrix &source_matrix)
	Copy constructor. More...

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

virtual	~CRMatrix ()
	Destructor, delete any allocated memory. More...

std::complex< double >	get_entry (const unsigned long &i, const unsigned long &j) const

std::complex< double > &	entry (const unsigned long &i, const unsigned long &j)
	The read-write access function is deliberately broken. More...

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

void	output_bottom_right_zero_helper (std::ostream &outfile) const

void	sparse_indexed_output_helper (std::ostream &outfile) const

void	clean_up_memory ()
	Wipe matrix data and set all values to 0. More...

void	build (const Vector< std::complex< double > > &value, const Vector< int > &column_index, const Vector< int > &row_start, const unsigned long &n, const unsigned long &m)

void	build_without_copy (std::complex< double > value, int column_index, int *row_start, const unsigned long &nnz, const unsigned long &n, const unsigned long &m)

Public Member Functions inherited from oomph::SparseMatrix< T, MATRIX_TYPE >
	SparseMatrix ()
	Default constructor. More...

	SparseMatrix (const SparseMatrix &source_matrix)
	Copy constructor. More...

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

virtual	~SparseMatrix ()
	Destructor, delete the memory associated with the values. More...

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

const T *	value () const
	Access to C-style value array (const version) 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...

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

Public Member Functions inherited from oomph::Matrix< T, MATRIX_TYPE >
	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...

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

T &	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::ComplexMatrixBase
	ComplexMatrixBase ()
	(Empty) constructor. More...

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

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

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

virtual void	solve (Vector< std::complex< double >> &rhs)

virtual void	solve (const Vector< std::complex< double >> &rhs, Vector< std::complex< double >> &soln)

virtual double	max_residual (const Vector< std::complex< double >> &x, const Vector< std::complex< double >> &rhs)

Protected Attributes
bool	Doc_stats_during_solve

Protected Attributes inherited from oomph::CRMatrix< std::complex< double > >
int *	Column_index
	Column index. More...

int *	Row_start
	Start index for row. More...

Protected Attributes inherited from oomph::SparseMatrix< T, MATRIX_TYPE >
T *	Value
	Internal representation of the matrix values, a pointer. More...

unsigned long	N
	Number of rows. More...

unsigned long	M
	Number of columns. More...

unsigned long	Nnz
	Number of non-zero values (i.e. size of Value array) More...

Private Attributes
void *	F_factors
	Storage for the LU factors as required by SuperLU. More...

int	Info
	Info flag for the SuperLU solver. More...

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

Static Protected Attributes inherited from oomph::SparseMatrix< T, MATRIX_TYPE >
static T	Zero = T(0)
	Dummy zero. More...

Detailed Description

A class for compressed row matrices.

Constructor & Destructor Documentation

◆ CRComplexMatrix() [1/3]

oomph::CRComplexMatrix::CRComplexMatrix ( )

inline

Default constructor.

                       : CRMatrix<std::complex<double>>(), F_factors(0), Info(0)
     {
       // By default SuperLU solves linear systems quietly
       Doc_stats_during_solve = false;
     }

References Doc_stats_during_solve.

◆ CRComplexMatrix() [2/3]

oomph::CRComplexMatrix::CRComplexMatrix	(	const Vector< std::complex< double >> &	value,
		const Vector< int > &	column_index,
		const Vector< int > &	row_start,
		const unsigned long &	n,
		const unsigned long &	m
	)

inline

Constructor: Pass vector of values, vector of column indices, vector of row starts and number of columns (can be suppressed for square matrices)

       : CRMatrix<std::complex<double>>(value, column_index, row_start, n, m),
         F_factors(0),
         Info(0)
     {
       // By default SuperLU solves linear systems quietly
       Doc_stats_during_solve = false;
     }

References Doc_stats_during_solve.

◆ CRComplexMatrix() [3/3]

oomph::CRComplexMatrix::CRComplexMatrix ( const CRComplexMatrix & matrix )

delete

Broken copy constructor.

◆ ~CRComplexMatrix()

virtual oomph::CRComplexMatrix::~CRComplexMatrix ( )

inlinevirtual

Destructor: Kill the LU decomposition if it has been computed.

     {
       clean_up_memory();
     }

References clean_up_memory().

Member Function Documentation

◆ clean_up_memory()

void oomph::CRComplexMatrix::clean_up_memory ( )

LU clean up (perhaps this should happen in the destructor)

Cleanup memory.

   {
     // Only bother if we've done an LU solve
     if (F_factors != 0)
     {
       // SuperLU expects compressed column format: Set flag for
       // transpose (0/1) = (false/true)
       int transpose = 1;
  
       // Doc (0/1) = (false/true)
       int doc = 0;
       if (Doc_stats_during_solve) doc = 1;
  
       // Copies so that const-ness is maintained
       int n_aux = int(N);
       int nnz_aux = int(Nnz);
  
       // SuperLU in three steps (lu decompose, back subst, cleanup)
       // Flag to indicate which solve step to do (1, 2 or 3)
       int i = 3;
       superlu_complex(&i,
                       &n_aux,
                       &nnz_aux,
                       0,
                       Value,
                       Column_index,
                       Row_start,
                       0,
                       &n_aux,
                       &transpose,
                       &doc,
                       &F_factors,
                       &Info);
     }
   }

References oomph::CRMatrix< std::complex< double > >::Column_index, Doc_stats_during_solve, F_factors, i, Info, int(), oomph::SparseMatrix< T, MATRIX_TYPE >::N, oomph::SparseMatrix< T, MATRIX_TYPE >::Nnz, oomph::CRMatrix< std::complex< double > >::Row_start, oomph::superlu_complex(), anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), and oomph::SparseMatrix< T, MATRIX_TYPE >::Value.

Referenced by ~CRComplexMatrix().

◆ disable_doc_stats()

void oomph::CRComplexMatrix::disable_doc_stats ( )

inline

the solve

     {
       Doc_stats_during_solve = false;
     }

References Doc_stats_during_solve.

◆ enable_doc_stats()

void oomph::CRComplexMatrix::enable_doc_stats ( )

inline

Set flag to indicate that stats are to be displayed during solution of linear system with SuperLU

     {
       Doc_stats_during_solve = true;
     }

References Doc_stats_during_solve.

◆ lubksub()

void oomph::CRComplexMatrix::lubksub ( Vector< std::complex< double >> & rhs )

virtual

LU back solve for given RHS.

Do back-substitution.

RHS vector

Reimplemented from oomph::ComplexMatrixBase.

   {
 #ifdef PARANOID
     if (N != rhs.size())
     {
       std::ostringstream error_message_stream;
       error_message_stream << "Size of RHS doesn't match matrix size\n"
                            << "N, rhs.size() " << N << " " << rhs.size()
                            << std::endl;
  
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     if (N != M)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "Can only solve for quadratic matrices\n"
                            << "N, M " << N << " " << M << std::endl;
  
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     std::complex<double>* b = new std::complex<double>[N];
  
     // Copy across
     for (unsigned long i = 0; i < N; i++)
     {
       b[i] = rhs[i];
     }
  
     // SuperLU expects compressed column format: Set flag for
     // transpose (0/1) = (false/true)
     int transpose = 1;
  
     // Doc (0/1) = (false/true)
     int doc = 0;
     if (Doc_stats_during_solve) doc = 1;
  
     // Number of RHSs
     int nrhs = 1;
  
     // Copies so that const-ness is maintained
     int n_aux = int(N);
     int nnz_aux = int(Nnz);
  
     // SuperLU in three steps (lu decompose, back subst, cleanup)
     // Do the solve phase
     int i = 2;
     superlu_complex(&i,
                     &n_aux,
                     &nnz_aux,
                     &nrhs,
                     Value,
                     Column_index,
                     Row_start,
                     b,
                     &n_aux,
                     &transpose,
                     &doc,
                     &F_factors,
                     &Info);
  
     // Copy b into rhs vector
     for (unsigned long i = 0; i < N; i++)
     {
       rhs[i] = b[i];
     }
  
     // Cleanup
     delete[] b;
   }

References b, oomph::CRMatrix< std::complex< double > >::Column_index, Doc_stats_during_solve, F_factors, i, Info, int(), oomph::SparseMatrix< T, MATRIX_TYPE >::N, oomph::SparseMatrix< T, MATRIX_TYPE >::Nnz, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRMatrix< std::complex< double > >::Row_start, oomph::superlu_complex(), anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), and oomph::SparseMatrix< T, MATRIX_TYPE >::Value.

◆ ludecompose()

int oomph::CRComplexMatrix::ludecompose ( )

virtual

LU decomposition using SuperLU.

Do LU decomposition and return sign of determinant.

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

Reimplemented from oomph::ComplexMatrixBase.

   {
 #ifdef PARANOID
     if (N != M)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "Can only solve for quadratic matrices\n"
                            << "N, M " << N << " " << M << std::endl;
  
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // SuperLU expects compressed column format: Set flag for
     // transpose (0/1) = (false/true)
     int transpose = 1;
  
     // Doc (0/1) = (false/true)
     int doc = 0;
     if (Doc_stats_during_solve) doc = 1;
  
     // Copies so that const-ness is maintained
     int n_aux = int(N);
     int nnz_aux = int(Nnz);
  
     // Integer to hold the sign of the determinant
     int sign = 0;
  
     // Just do the lu decompose phase (i=1)
     int i = 1;
     sign = superlu_complex(&i,
                            &n_aux,
                            &nnz_aux,
                            0,
                            Value,
                            Column_index,
                            Row_start,
                            0,
                            &n_aux,
                            &transpose,
                            &doc,
                            &F_factors,
                            &Info);
     // Return sign of determinant
     return sign;
   }

References oomph::CRMatrix< std::complex< double > >::Column_index, Doc_stats_during_solve, F_factors, i, Info, int(), oomph::SparseMatrix< T, MATRIX_TYPE >::N, oomph::SparseMatrix< T, MATRIX_TYPE >::Nnz, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRMatrix< std::complex< double > >::Row_start, SYCL::sign(), oomph::superlu_complex(), anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), and oomph::SparseMatrix< T, MATRIX_TYPE >::Value.

◆ multiply()

void oomph::CRComplexMatrix::multiply	(	const Vector< std::complex< double >> &	x,
		Vector< std::complex< double >> &	soln
	)

virtual

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

Multiply the matrix by the vector x.

Implements oomph::ComplexMatrixBase.

   {
 #ifdef PARANOID
     // Check to see x.size() = ncol()
     if (x.size() != M)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The x vector is not the right size. It is "
                            << x.size() << ", it should be " << M << std::endl;
  
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     if (soln.size() != N)
     {
       // Resize and initialize the solution vector
       soln.resize(N);
     }
     for (unsigned long i = 0; i < N; i++)
     {
       soln[i] = 0.0;
       for (long k = Row_start[i]; k < Row_start[i + 1]; k++)
       {
         unsigned long j = Column_index[k];
         std::complex<double> a_ij = Value[k];
         soln[i] += a_ij * x[j];
       }
     }
   }

References oomph::CRMatrix< std::complex< double > >::Column_index, i, j, k, oomph::SparseMatrix< T, MATRIX_TYPE >::N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRMatrix< std::complex< double > >::Row_start, oomph::SparseMatrix< T, MATRIX_TYPE >::Value, and plotDoE::x.

◆ multiply_transpose()

void oomph::CRComplexMatrix::multiply_transpose	(	const Vector< std::complex< double >> &	x,
		Vector< std::complex< double >> &	soln
	)

virtual

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

Implements oomph::ComplexMatrixBase.

   {
 #ifdef PARANOID
     // Check to see x.size() = nrow()
     if (x.size() != N)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The x vector is not the right size. It is "
                            << x.size() << ", it should be " << N << std::endl;
  
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     if (soln.size() != M)
     {
       // Resize and initialize the solution vector
       soln.resize(M);
     }
  
     // Initialise the solution
     for (unsigned long i = 0; i < M; i++)
     {
       soln[i] = 0.0;
     }
  
     // 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];
         std::complex<double> a_ij = Value[k];
         soln[j] += a_ij * x[i];
       }
     }
   }

References oomph::CRMatrix< std::complex< double > >::Column_index, i, j, k, oomph::SparseMatrix< T, MATRIX_TYPE >::M, oomph::SparseMatrix< T, MATRIX_TYPE >::N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRMatrix< std::complex< double > >::Row_start, oomph::SparseMatrix< T, MATRIX_TYPE >::Value, and plotDoE::x.

◆ ncol()

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

inlinevirtual

Return the number of columns of the matrix.

Implements oomph::ComplexMatrixBase.

     {
       return CRMatrix<std::complex<double>>::ncol();
     }

◆ nrow()

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

inlinevirtual

Return the number of rows of the matrix.

Implements oomph::ComplexMatrixBase.

     {
       return CRMatrix<std::complex<double>>::nrow();
     }

◆ operator()()

std::complex<double> oomph::CRComplexMatrix::operator()	(	const unsigned long &	i,
		const unsigned long &	j
	)		const

inlinevirtual

Overload the round-bracket access operator for read-only access.

Implements oomph::ComplexMatrixBase.

     {
       return CRMatrix<std::complex<double>>::get_entry(i, j);
     }

References oomph::CRMatrix< std::complex< double > >::get_entry(), i, and j.

◆ operator=()

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

delete

Broken assignment operator.

◆ residual()

void oomph::CRComplexMatrix::residual	(	const Vector< std::complex< double >> &	x,
		const Vector< std::complex< double >> &	b,
		Vector< std::complex< double >> &	residual
	)

virtual

Find the residual to x of Ax=b, i.e. r=b-Ax.

Find the residulal to x of Ax=b, ie r=b-Ax.

Implements oomph::ComplexMatrixBase.

   {
 #ifdef PARANOID
     // Check that size of rhs = nrow()
     if (rhs.size() != N)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The rhs vector is not the right size. It is "
                            << rhs.size() << ", it should be " << N << std::endl;
  
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
     // Check that the size of x is correct
     if (x.size() != M)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "The x vector is not the right size. It is "
                            << x.size() << ", it should be " << M << std::endl;
  
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     if (residual.size() != N)
     {
       residual.resize(N);
     }
  
     for (unsigned long i = 0; i < N; i++)
     {
       residual[i] = rhs[i];
       for (long k = Row_start[i]; k < Row_start[i + 1]; k++)
       {
         unsigned long j = Column_index[k];
         std::complex<double> a_ij = Value[k];
         residual[i] -= a_ij * x[j];
       }
     }
   }

References oomph::CRMatrix< std::complex< double > >::Column_index, i, j, k, oomph::SparseMatrix< T, MATRIX_TYPE >::N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRMatrix< std::complex< double > >::Row_start, oomph::SparseMatrix< T, MATRIX_TYPE >::Value, and plotDoE::x.

Member Data Documentation

◆ Doc_stats_during_solve

bool oomph::CRComplexMatrix::Doc_stats_during_solve

protected

Flag to indicate if stats are to be displayed during solution of linear system with SuperLU

Referenced by clean_up_memory(), CRComplexMatrix(), disable_doc_stats(), enable_doc_stats(), lubksub(), and ludecompose().

◆ F_factors

void* oomph::CRComplexMatrix::F_factors

private

Storage for the LU factors as required by SuperLU.

Referenced by clean_up_memory(), lubksub(), and ludecompose().

◆ Info

int oomph::CRComplexMatrix::Info

private

Info flag for the SuperLU solver.

Referenced by clean_up_memory(), lubksub(), and ludecompose().

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

Public Member Functions

Protected Attributes

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ CRComplexMatrix() [1/3]

◆ CRComplexMatrix() [2/3]

◆ CRComplexMatrix() [3/3]

◆ ~CRComplexMatrix()

Member Function Documentation

◆ clean_up_memory()

◆ disable_doc_stats()

◆ enable_doc_stats()

◆ lubksub()

◆ ludecompose()

◆ multiply()

◆ multiply_transpose()

◆ ncol()

◆ nrow()

◆ operator()()

◆ operator=()

◆ residual()

Member Data Documentation

◆ Doc_stats_during_solve

◆ F_factors

◆ Info