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

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

RSS