oomph::CCDoubleMatrix Class Reference

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

#include <matrices.h>

Inheritance diagram for oomph::CCDoubleMatrix:

Public Member Functions
	CCDoubleMatrix ()
	Default constructor. More...
	CCDoubleMatrix (const Vector< double > &value, const Vector< int > &row_index_, const Vector< int > &column_start_, const unsigned long &n, const unsigned long &m)
	CCDoubleMatrix (const CCDoubleMatrix &matrix)=delete
	Broken copy constructor. More...
void	operator= (const CCDoubleMatrix &)=delete
	Broken assignment operator. More...
virtual	~CCDoubleMatrix ()
	Destructor: Kill the LU factors if they have been setup. 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...
double	operator() (const unsigned long &i, const unsigned long &j) const
virtual void	ludecompose ()
	LU decomposition using SuperLU. 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 CCDoubleMatrix &matrix_in, CCDoubleMatrix &result)
void	matrix_reduction (const double &alpha, CCDoubleMatrix &reduced_matrix)
unsigned &	matrix_matrix_multiply_method ()
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::CCMatrix< double >
	CCMatrix ()
	Default constructor. More...
	CCMatrix (const Vector< 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...
double	get_entry (const unsigned long &i, const unsigned long &j) const
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< double > &value, const Vector< int > &row_index, const Vector< int > &column_start, const unsigned long &n, const unsigned long &m)
void	build_without_copy (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	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

Private Attributes
unsigned	Matrix_matrix_multiply_method
	Flag to determine which matrix-matrix multiplication method is used. 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
Protected Attributes inherited from oomph::DoubleMatrixBase
LinearSolver *	Linear_solver_pt
LinearSolver *	Default_linear_solver_pt
Protected Attributes inherited from oomph::CCMatrix< 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...
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

CCDoubleMatrix() [1/3]

oomph::CCDoubleMatrix::CCDoubleMatrix

(

)

Default constructor.

//////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// Default constructor, set the default linear solver and matrix-matrix multiplication method.

                                 : CCMatrix<double>()
  {
    Linear_solver_pt = Default_linear_solver_pt = new SuperLUSolver;
    Matrix_matrix_multiply_method = 2;
  }

References oomph::DoubleMatrixBase::Default_linear_solver_pt, oomph::DoubleMatrixBase::Linear_solver_pt, and Matrix_matrix_multiply_method.

CCDoubleMatrix() [2/3]

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

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<double>(value, row_index, column_start, n, m)
  {
    Linear_solver_pt = Default_linear_solver_pt = new SuperLUSolver;
    Matrix_matrix_multiply_method = 2;
  }

References oomph::DoubleMatrixBase::Default_linear_solver_pt, oomph::DoubleMatrixBase::Linear_solver_pt, and Matrix_matrix_multiply_method.

CCDoubleMatrix() [3/3]

oomph::CCDoubleMatrix::CCDoubleMatrix ( const CCDoubleMatrix &  matrix )

delete

Broken copy constructor.

~CCDoubleMatrix()

oomph::CCDoubleMatrix::~CCDoubleMatrix ( )

virtual

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

Destructor: delete the default linear solver.

  {
    delete Default_linear_solver_pt;
  }

References oomph::DoubleMatrixBase::Default_linear_solver_pt.

Member Function Documentation

lubksub()

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

virtual

LU back solve for given RHS.

Do the backsubstitution.

  {
    static_cast<SuperLUSolver*>(Default_linear_solver_pt)->backsub(rhs, rhs);
  }

References oomph::DoubleMatrixBase::Default_linear_solver_pt.

ludecompose()

void oomph::CCDoubleMatrix::ludecompose ( )

virtual

LU decomposition using SuperLU.

Perform LU decomposition. Return the sign of the determinant.

  {
    static_cast<SuperLUSolver*>(Default_linear_solver_pt)->factorise(this);
  }

References oomph::DoubleMatrixBase::Default_linear_solver_pt.

matrix_matrix_multiply_method()

unsigned& oomph::CCDoubleMatrix::matrix_matrix_multiply_method ( )

inline

Access function to Matrix_matrix_multiply_method, the flag which determines the matrix matrix multiplication method used. 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...

    {
      return Matrix_matrix_multiply_method;
    }

References Matrix_matrix_multiply_method.

Referenced by main().

matrix_reduction()

void oomph::CCDoubleMatrix::matrix_reduction	(	const double &	alpha,
		CCDoubleMatrix &	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 columns in matrix
    long n_coln = ncol();
 
    Vector<double> max_row(nrow(), 0.0);
 
    // Here's the packed format for the new matrix
    Vector<int> B_row_start(1);
    Vector<int> B_column_index;
    Vector<double> B_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 columns
    for (long i = 0; i < n_coln; i++)
    {
      // Loop over entries in columns
      for (long j = Column_start[i]; j < Column_start[i + 1]; j++)
      {
        // Find max. value in row
        if (std::fabs(Value[j]) > max_row[Row_index[j]])
        {
          max_row[Row_index[j]] = std::fabs(Value[j]);
        }
      }
 
      // Decide if we need to retain the entries in the row
      for (long j = Column_start[i]; j < Column_start[i + 1]; j++)
      {
        // If we're on the diagonal or the value is sufficiently large: retain
        // i.e. copy across.
        if (i == Row_index[j] ||
            std::fabs(Value[j]) > alpha * max_row[Row_index[j]])
        {
          B_value.push_back(Value[j]);
          B_column_index.push_back(Row_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<CCDoubleMatrix&>(reduced_matrix)
      .build(B_value, B_column_index, B_row_start, nrow(), ncol());
  }

References alpha, oomph::CCMatrix< double >::build(), oomph::CCMatrix< double >::Column_start, boost::multiprecision::fabs(), i, j, k, ncol(), nrow(), oomph::CCMatrix< double >::Row_index, and oomph::SparseMatrix< T, MATRIX_TYPE >::Value.

multiply() [1/2]

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

Function to multiply this matrix by the CCDoubleMatrix matrix_in The multiplication method used can be selected using the flag Matrix_matrix_multiply_method. By default Method 2 is used. 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...

  {
#ifdef PARANOID
    // check matrix dimensions are compatible
    if (this->ncol() != matrix_in.nrow())
    {
      std::ostringstream error_message;
      error_message
        << "Matrix dimensions incompatable for matrix-matrix multiplication"
        << "ncol() for first matrix:" << this->ncol()
        << "nrow() for second matrix: " << matrix_in.nrow();
 
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // 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* Column_start;
    double* Value;
    int* Row_index;
 
    // get pointers to matrix_in
    const int* matrix_in_col_start = matrix_in.column_start();
    const int* matrix_in_row_index = matrix_in.row_index();
    const double* matrix_in_value = matrix_in.value();
 
    // get pointers to this matrix
    const double* this_value = this->value();
    const int* this_col_start = this->column_start();
    const int* this_row_index = this->row_index();
 
    // set method
    unsigned method = Matrix_matrix_multiply_method;
 
    // clock_t clock1 = clock();
 
    // METHOD 1
    // --------
    if (method == 1)
    {
      // allocate storage for column starts
      Column_start = new int[M + 1];
      Column_start[0] = 0;
 
      // a set to store number of non-zero rows in each column of result
      std::set<unsigned> rows;
 
      // run through columns of this matrix and matrix_in to find number of
      // non-zero entries in each column of result
      for (unsigned long this_col = 0; this_col < M; this_col++)
      {
        // run through non-zeros in this_col of this matrix
        for (int this_ptr = this_col_start[this_col];
             this_ptr < this_col_start[this_col + 1];
             this_ptr++)
        {
          // find row index for non-zero
          unsigned matrix_in_col = this_row_index[this_ptr];
 
          // run through corresponding column in matrix_in
          for (int matrix_in_ptr = matrix_in_col_start[matrix_in_col];
               matrix_in_ptr < matrix_in_col_start[matrix_in_col + 1];
               matrix_in_ptr++)
          {
            // find row index for non-zero in matrix_in and store in rows
            rows.insert(matrix_in_row_index[matrix_in_ptr]);
          }
        }
        // update Column_start
        Column_start[this_col + 1] = Column_start[this_col] + rows.size();
 
        // wipe values in rows
        rows.clear();
      }
 
      // set Nnz
      Nnz = Column_start[M];
 
      // allocate arrays for result
      Value = new double[Nnz];
      Row_index = new int[Nnz];
 
      // set all values of Row_index to -1
      for (unsigned long i = 0; i < Nnz; i++) Row_index[i] = -1;
 
      // Calculate values for result - first run through columns of this matrix
      for (unsigned long this_col = 0; this_col < M; this_col++)
      {
        // run through non-zeros in this_column
        for (int this_ptr = this_col_start[this_col];
             this_ptr < this_col_start[this_col + 1];
             this_ptr++)
        {
          // find value of non-zero
          double this_val = this_value[this_ptr];
 
          // find row associated with non-zero
          unsigned matrix_in_col = this_row_index[this_ptr];
 
          // run through corresponding column in matrix_in
          for (int matrix_in_ptr = matrix_in_col_start[matrix_in_col];
               matrix_in_ptr < matrix_in_col_start[matrix_in_col + 1];
               matrix_in_ptr++)
          {
            // find row index for non-zero in matrix_in
            int row = matrix_in_row_index[matrix_in_ptr];
 
            // find position in result to insert value
            for (int ptr = Column_start[this_col];
                 ptr <= Column_start[this_col + 1];
                 ptr++)
            {
              if (ptr == Column_start[this_col + 1])
              {
                // error - have passed end of column without finding
                // correct row index
                std::ostringstream error_message;
                error_message << "Error inserting value in result";
 
                throw OomphLibError(error_message.str(),
                                    OOMPH_CURRENT_FUNCTION,
                                    OOMPH_EXCEPTION_LOCATION);
              }
              else if (Row_index[ptr] == -1)
              {
                // first entry for this row index
                Row_index[ptr] = row;
                Value[ptr] = this_val * matrix_in_value[matrix_in_ptr];
                break;
              }
              else if (Row_index[ptr] == row)
              {
                // row 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>[M];
 
      // run through columns of this matrix
      for (unsigned long this_col = 0; this_col < M; this_col++)
      {
        // run through non-zeros in this_col
        for (int this_ptr = this_col_start[this_col];
             this_ptr < this_col_start[this_col + 1];
             this_ptr++)
        {
          // find value of non-zero
          double this_val = this_value[this_ptr];
 
          // find row index associated with non-zero
          unsigned matrix_in_col = this_row_index[this_ptr];
 
          // run through corresponding column in matrix_in
          for (int matrix_in_ptr = matrix_in_col_start[matrix_in_col];
               matrix_in_ptr < matrix_in_col_start[matrix_in_col + 1];
               matrix_in_ptr++)
          {
            // find row index for non-zero in matrix_in
            int row = matrix_in_row_index[matrix_in_ptr];
 
            // insert value
            result_maps[this_col][row] +=
              this_val * matrix_in_value[matrix_in_ptr];
          }
        }
      }
 
      // allocate Column_start
      Column_start = new int[M + 1];
 
      // copy across column starts
      Column_start[0] = 0;
      for (unsigned long col = 0; col < M; col++)
      {
        int size = result_maps[col].size();
        Column_start[col + 1] = Column_start[col] + size;
      }
 
      // set Nnz
      Nnz = Column_start[M];
 
      // allocate other arrays
      Value = new double[Nnz];
      Row_index = new int[Nnz];
 
      // copy values and row indices
      for (unsigned long col = 0; col < M; col++)
      {
        unsigned ptr = Column_start[col];
        for (std::map<int, double>::iterator i = result_maps[col].begin();
             i != result_maps[col].end();
             i++)
        {
          Row_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_rows(N);
      std::vector<std::vector<double>> result_vals(N);
 
      // run through the columns of this matrix
      for (unsigned long this_col = 0; this_col < M; this_col++)
      {
        // run through non-zeros in this_col
        for (int this_ptr = this_col_start[this_col];
             this_ptr < this_col_start[this_col + 1];
             this_ptr++)
        {
          // find value of non-zero
          double this_val = this_value[this_ptr];
 
          // find row index associated with non-zero
          unsigned matrix_in_col = this_row_index[this_ptr];
 
          // run through corresponding column in matrix_in
          for (int matrix_in_ptr = matrix_in_col_start[matrix_in_col];
               matrix_in_ptr < matrix_in_col_start[matrix_in_col + 1];
               matrix_in_ptr++)
          {
            // find row index for non-zero in matrix_in
            int row = matrix_in_row_index[matrix_in_ptr];
 
            // insert value
            int size = result_rows[this_col].size();
            for (int i = 0; i <= size; i++)
            {
              if (i == size)
              {
                // first entry for this row index
                result_rows[this_col].push_back(row);
                result_vals[this_col].push_back(this_val *
                                                matrix_in_value[matrix_in_ptr]);
              }
              else if (row == result_rows[this_col][i])
              {
                // row index already exists
                result_vals[this_col][i] +=
                  this_val * matrix_in_value[matrix_in_ptr];
                break;
              }
            }
          }
        }
      }
 
      // allocate Column_start
      Column_start = new int[M + 1];
 
      // copy across column starts
      Column_start[0] = 0;
      for (unsigned long col = 0; col < M; col++)
      {
        int size = result_rows[col].size();
        Column_start[col + 1] = Column_start[col] + size;
      }
 
      // set Nnz
      Nnz = Column_start[M];
 
      // allocate other arrays
      Value = new double[Nnz];
      Row_index = new int[Nnz];
 
      // copy across values and row indices
      for (unsigned long col = 0; col < N; col++)
      {
        unsigned ptr = Column_start[col];
        unsigned n_rows = result_rows[col].size();
        for (unsigned i = 0; i < n_rows; i++)
        {
          Row_index[ptr] = result_rows[col][i];
          Value[ptr] = result_vals[col][i];
          ptr++;
        }
      }
    }
 
    // INCORRECT VALUE FOR METHOD
    else
    {
      std::ostringstream error_message;
      error_message << "Incorrect method set in matrix-matrix multiply"
                    << "method=" << method << " not allowed";
 
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    result.build_without_copy(Value, Row_index, Column_start, Nnz, N, M);
  }

References oomph::CCMatrix< T >::build_without_copy(), col(), oomph::CCMatrix< double >::column_start(), oomph::CCMatrix< T >::column_start(), oomph::CCMatrix< double >::Column_start, i, oomph::SparseMatrix< T, MATRIX_TYPE >::M, Matrix_matrix_multiply_method, oomph::SparseMatrix< T, MATRIX_TYPE >::N, ncol(), oomph::SparseMatrix< T, MATRIX_TYPE >::Nnz, nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, row(), oomph::CCMatrix< T >::row_index(), oomph::CCMatrix< double >::row_index(), oomph::CCMatrix< double >::Row_index, rows, size, oomph::SparseMatrix< T, MATRIX_TYPE >::Value, and oomph::SparseMatrix< T, MATRIX_TYPE >::value().

multiply() [2/2]

void oomph::CCDoubleMatrix::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 to see if x.size() = ncol().
    if (x.nrow() != this->ncol())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The x vector is not the right size. It is "
                           << x.nrow() << ", it should be " << this->ncol()
                           << std::endl;
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // check that x is not distributed
    if (x.distributed())
    {
      std::ostringstream error_message_stream;
      error_message_stream
        << "The x vector cannot be distributed for CCDoubleMatrix "
        << "matrix-vector multiply" << std::endl;
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // if soln is setup...
    if (soln.built())
    {
      // check that soln is not distributed
      if (soln.distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The x vector cannot be distributed for CCDoubleMatrix "
          << "matrix-vector multiply" << std::endl;
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      if (soln.nrow() != this->nrow())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The soln vector is setup and therefore must have the same "
          << "number of rows as the matrix";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      if (*soln.distribution_pt()->communicator_pt() !=
          *x.distribution_pt()->communicator_pt())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The soln vector and the x vector must have the same communicator"
          << std::endl;
        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->nrow(), false);
      soln.build(dist_pt, 0.0);
      delete dist_pt;
    }
 
    // zero
    soln.initialise(0.0);
 
    // multiply
    double* soln_pt = soln.values_pt();
    const double* x_pt = x.values_pt();
    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];
        double a_ij = Value[k];
        soln_pt[i] += a_ij * x_pt[j];
      }
    }
  }

References oomph::DoubleVector::build(), oomph::DoubleVector::built(), oomph::CCMatrix< double >::Column_start, oomph::LinearAlgebraDistribution::communicator_pt(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), i, oomph::DoubleVector::initialise(), j, k, oomph::SparseMatrix< T, MATRIX_TYPE >::N, ncol(), oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CCMatrix< double >::Row_index, oomph::SparseMatrix< T, MATRIX_TYPE >::Value, oomph::DoubleVector::values_pt(), and plotDoE::x.

Referenced by main().

multiply_transpose()

void oomph::CCDoubleMatrix::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 to see if x.size() = ncol().
    if (x.nrow() != this->nrow())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The x vector is not the right size. It is "
                           << x.nrow() << ", it should be " << this->nrow()
                           << std::endl;
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // check that x is not distributed
    if (x.distributed())
    {
      std::ostringstream error_message_stream;
      error_message_stream
        << "The x vector cannot be distributed for CCDoubleMatrix "
        << "matrix-vector multiply" << std::endl;
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // if soln is setup...
    if (soln.built())
    {
      // check that soln is not distributed
      if (soln.distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The x vector cannot be distributed for CCDoubleMatrix "
          << "matrix-vector multiply" << std::endl;
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      if (soln.nrow() != this->ncol())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The soln vector is setup and therefore must have the same "
          << "number of columns as the matrix";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      if (*soln.distribution_pt()->communicator_pt() !=
          *x.distribution_pt()->communicator_pt())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The soln vector and the x vector must have the same communicator"
          << std::endl;
        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(), false);
      soln.build(dist_pt, 0.0);
      delete dist_pt;
    }
 
    // zero
    soln.initialise(0.0);
 
    // Matrix vector product
    double* soln_pt = soln.values_pt();
    const double* x_pt = x.values_pt();
    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];
        double a_ij = Value[k];
        soln_pt[j] += a_ij * x_pt[i];
      }
    }
  }

References oomph::DoubleVector::build(), oomph::DoubleVector::built(), oomph::CCMatrix< double >::Column_start, oomph::LinearAlgebraDistribution::communicator_pt(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), i, oomph::DoubleVector::initialise(), j, k, oomph::SparseMatrix< T, MATRIX_TYPE >::N, oomph::DistributableLinearAlgebraObject::nrow(), nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CCMatrix< double >::Row_index, oomph::SparseMatrix< T, MATRIX_TYPE >::Value, oomph::DoubleVector::values_pt(), and plotDoE::x.

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

ncol()

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

inlinevirtual

Return the number of columns of the matrix.

Implements oomph::DoubleMatrixBase.

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

References oomph::SparseMatrix< T, CCMatrix< T > >::ncol().

Referenced by matrix_reduction(), and multiply().

nrow()

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

inlinevirtual

Return the number of rows of the matrix.

Implements oomph::DoubleMatrixBase.

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

References oomph::SparseMatrix< T, CCMatrix< T > >::nrow().

Referenced by matrix_reduction(), multiply(), multiply_transpose(), and oomph::ILUZeroPreconditioner< CCDoubleMatrix >::setup().

operator()()

double oomph::CCDoubleMatrix::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::DoubleMatrixBase.

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

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

operator=()

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

delete

Broken assignment operator.

Member Data Documentation

Matrix_matrix_multiply_method

unsigned oomph::CCDoubleMatrix::Matrix_matrix_multiply_method

private

Flag to determine which matrix-matrix multiplication method is used.

Referenced by CCDoubleMatrix(), matrix_matrix_multiply_method(), and multiply().

---

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

• matrices.h
• matrices.cc