A class for compressed column matrices that store doubles. More...

#include <complex_matrices.h>

Inheritance diagram for oomph::CCComplexMatrix:

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

	CCComplexMatrix (const Vector< std::complex< double >> &value, const Vector< int > &row_index, const Vector< int > &column_start, const unsigned long &n, const unsigned long &m)

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

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

virtual	~CCComplexMatrix ()
	Destructor: Kill the LU factors if they have been setup. 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

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 residulal to x of Ax=b, ie 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::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)

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

	CCMatrix (const Vector< std::complex< double > > &value, const Vector< int > &row_index_, const Vector< int > &column_start_, const unsigned long &n, const unsigned long &m)

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

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

virtual	~CCMatrix ()
	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)

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

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

int *	row_index ()
	Access to C-style row index array. More...

const int *	row_index () const
	Access to C-style row 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 > &row_index, const Vector< int > &column_start, const unsigned long &n, const unsigned long &m)

void	build_without_copy (std::complex< double > value, int row_index, int *column_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

Protected Attributes
bool	Doc_stats_during_solve

Protected Attributes inherited from oomph::CCMatrix< std::complex< double > >
int *	Row_index
	Row index. More...

int *	Column_start
	Start index for column. 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 column matrices that store doubles.

Constructor & Destructor Documentation

◆ CCComplexMatrix() [1/3]

oomph::CCComplexMatrix::CCComplexMatrix ( )

inline

Default constructor.

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

◆ CCComplexMatrix() [2/3]

oomph::CCComplexMatrix::CCComplexMatrix	(	const Vector< std::complex< double >> &	value,
		const Vector< int > &	row_index,
		const Vector< int > &	column_start,
		const unsigned long &	n,
		const unsigned long &	m
	)

inline

Constructor: Pass vector of values, vector of row indices, vector of column starts and number of rows (can be suppressed for square matrices). Number of nonzero entries is read off from value, so make sure the vector has been shrunk to its correct length.

       : CCMatrix<std::complex<double>>(value, row_index, column_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.

◆ CCComplexMatrix() [3/3]

oomph::CCComplexMatrix::CCComplexMatrix ( const CCComplexMatrix & matrix )

delete

Broken copy constructor.

◆ ~CCComplexMatrix()

virtual oomph::CCComplexMatrix::~CCComplexMatrix ( )

inlinevirtual

Destructor: Kill the LU factors if they have been setup.

     {
       clean_up_memory();
     }

References clean_up_memory().

Member Function Documentation

◆ clean_up_memory()

void oomph::CCComplexMatrix::clean_up_memory ( )

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

Cleanup storage.

   {
     // 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 = 0;
  
       // Doc (0/1) = (false/true)
       int doc = 0;
       if (Doc_stats_during_solve) doc = 1;
  
  
       // Cast to integers for stupid SuperLU
       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,
                       Row_index,
                       Column_start,
                       0,
                       &n_aux,
                       &transpose,
                       &doc,
                       &F_factors,
                       &Info);
     }
   }

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

Referenced by ~CCComplexMatrix().

◆ disable_doc_stats()

void oomph::CCComplexMatrix::disable_doc_stats ( )

inline

the solve

     {
       Doc_stats_during_solve = false;
     }

References Doc_stats_during_solve.

◆ enable_doc_stats()

void oomph::CCComplexMatrix::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::CCComplexMatrix::lubksub ( Vector< std::complex< double >> & rhs )

virtual

LU back solve for given RHS.

Do the backsubstitution.

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 = 0;
  
     // Doc (0/1) = (false/true)
     int doc = 0;
     if (Doc_stats_during_solve) doc = 1;
  
     // Number of RHSs
     int nrhs = 1;
  
  
     // Cast to integers for stupid SuperLU
     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,
                     Row_index,
                     Column_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::CCMatrix< std::complex< double > >::Column_start, 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::CCMatrix< std::complex< double > >::Row_index, oomph::superlu_complex(), anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), and oomph::SparseMatrix< T, MATRIX_TYPE >::Value.

◆ ludecompose()

int oomph::CCComplexMatrix::ludecompose ( )

virtual

LU decomposition using SuperLU.

Perform LU decomposition. Return the sign of the 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 = 0;
  
     // Doc (0/1) = (false/true)
     int doc = 0;
     if (Doc_stats_during_solve) doc = 1;
  
     // Cast to integers for stupid SuperLU
     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,
                            Row_index,
                            Column_start,
                            0,
                            &n_aux,
                            &transpose,
                            &doc,
                            &F_factors,
                            &Info);
  
     // Return the sign of the determinant
     return sign;
   }

References oomph::CCMatrix< std::complex< double > >::Column_start, 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::CCMatrix< std::complex< double > >::Row_index, SYCL::sign(), oomph::superlu_complex(), anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), and oomph::SparseMatrix< T, MATRIX_TYPE >::Value.

◆ multiply()

void oomph::CCComplexMatrix::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 if 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 i = 0; i < N; i++)
     {
       soln[i] = 0.0;
     }
  
     for (unsigned long j = 0; j < N; j++)
     {
       for (long k = Column_start[j]; k < Column_start[j + 1]; k++)
       {
         unsigned long i = Row_index[k];
         std::complex<double> a_ij = Value[k];
         soln[i] += a_ij * x[j];
       }
     }
   }

References oomph::CCMatrix< std::complex< double > >::Column_start, i, j, k, oomph::SparseMatrix< T, MATRIX_TYPE >::N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CCMatrix< std::complex< double > >::Row_index, oomph::SparseMatrix< T, MATRIX_TYPE >::Value, and plotDoE::x.

◆ multiply_transpose()

void oomph::CCComplexMatrix::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 = Column_start[i]; k < Column_start[i + 1]; k++)
       {
         unsigned long j = Row_index[k];
         std::complex<double> a_ij = Value[k];
         soln[j] += a_ij * x[i];
       }
     }
   }

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

◆ ncol()

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

inlinevirtual

Return the number of columns of the matrix.

Implements oomph::ComplexMatrixBase.

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

◆ nrow()

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

inlinevirtual

Return the number of rows of the matrix.

Implements oomph::ComplexMatrixBase.

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

◆ operator()()

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

inlinevirtual

Overload the round-bracket access operator to provide read-only (const) access to the data

Implements oomph::ComplexMatrixBase.

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

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

◆ operator=()

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

delete

Broken assignment operator.

◆ residual()

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

virtual

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

Work out residual vector r = b-Ax for candidate solution x.

Implements oomph::ComplexMatrixBase.

   {
 #ifdef PARANOID
     // Check Matrix is square
     if (N != M)
     {
       throw OomphLibError(
         "This matrix is not square, the matrix MUST be square!",
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
     }
     // 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() != 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
  
     unsigned long r_n = residual.size();
     if (r_n != N)
     {
       residual.resize(N);
     }
  
     // Need to do this in loop over rows
     for (unsigned i = 0; i < N; i++)
     {
       residual[i] = rhs[i];
     }
     // Now loop over columns
     for (unsigned long j = 0; j < N; j++)
     {
       for (long k = Column_start[j]; k < Column_start[j + 1]; k++)
       {
         unsigned long i = Row_index[k];
         std::complex<double> a_ij = Value[k];
         residual[i] -= a_ij * x[j];
       }
     }
   }

References oomph::CCMatrix< std::complex< double > >::Column_start, i, j, k, oomph::SparseMatrix< T, MATRIX_TYPE >::N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CCMatrix< std::complex< double > >::Row_index, oomph::SparseMatrix< T, MATRIX_TYPE >::Value, and plotDoE::x.

Member Data Documentation

◆ Doc_stats_during_solve

bool oomph::CCComplexMatrix::Doc_stats_during_solve

protected

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

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

◆ F_factors

void* oomph::CCComplexMatrix::F_factors

private

Storage for the LU factors as required by SuperLU.

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

◆ Info

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

◆ CCComplexMatrix() [1/3]

◆ CCComplexMatrix() [2/3]

◆ CCComplexMatrix() [3/3]

◆ ~CCComplexMatrix()

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