The right-hand-side is usually dense or at least it was the case in my applications. In this case I have used the class Matrix that was already employed previously in sections devoted to dense matrices. I will describe this class in more detail here.

The class is declared in matrices.h and actually most functions are defined there as well. In C++ a one-dimensional array of doubles is well represented by vector<double> and the class Matrix is derived from it:

class Matrix : public vector<double>
{
  size_t m;
  size_t n;
public:

 The class has two private members that define the matrix dimension A(m, n) that could be accessed by the next two functions

  size_t NRows() const {return m;}
  size_t NCols() const {return n;}

The access function that overloads the parenthesis operator, A(i, j)

  double& operator()(size_t i, size_t j)
    {return operator[](i + j*m);}

defines that the matrix storage is assumed to be column-wise. In other world, in the underlying one-dimensional array first stores the first column of the matrix, then the second and so on.  This is accepted in Fortran and hence facilitates the use of Fortran subroutines. Note that the access operator from vector<double> can be also employed, for example A[i]. This could be useful when one is working with a vector A(m, 1) or when it is necessary to iterate through all the members of the matrix.

Another access operator

  double* column(size_t j)
    {return &*(begin() + j*m);}

returns the pointer to the i-th columns (i-th vector in the matrix). This could be sometime useful.

The functions

  size_t nnz()
    {return size() - count(begin(), end(), 0.);}
  size_t nnz(size_t j)
    {return NRows() - count(column(j), column(j) + NRows(), 0.);}

returns the number of nonzero entries in the matrix and in the i-th column accordingly.

The class has three constructors

  Matrix() : m(0), n(0) {}
  Matrix(size_t m_, size_t n_) :
    vector<double>(m_*n_), m(m_), n(n_) {}
  Matrix(size_t m_, size_t n_, double val) :
    vector<double>(m_*n_, val), m(m_), n(n_) {}

to create an empty matrix

Matrix A;

a matrix with a given dimension

Matrix A(m, n);

and a matrix with a given dimension and given value for all entries

Matrix A(m, n, 1);

The matrix could be dynamically resized A.resize(m1, n1)

  void resize(size_t m_, size_t n_)
    {m = m_; n = n_; vector<double>::resize(m_*n_);}

but if this is required often, then it may be more efficient first reserve the maximum matrix size, A.reserve(m_max, n_max), the product of m_max*n_max being the most important

  void reserve(size_t m_, size_t n_)
    {vector<double>::reserve(m_*n_);}

The function clear removes the content of the matrix, A.clear()

  void clear() {m = n = 0; vector<double>::clear();}

but it does not free memory back. If the matrix was huge and you need the memory for other purposes, please use A.clearMemory() instead.

  void clearMemory() {Matrix empty; swap(empty);}

Finally functions

  void read(istream &in);  
  void write(ostream &out);

allows us to read and write the matrix in the Matrix Market format. Their definitions are in matrices.cpp.

Previous

Introduction

Next

Class SparseMatrix
Sample code in C++ to run TAUCS
Sample code in C++ to run MUMPS
Class LinearSolver
  Class UMFPACK
  Class TAUCS
  Class MUMPS and PARDISO
  Gluing all together: LinearSolver::create
Sample driver to run a sparse solver