When one programs a dense matrix, basically there is necessary only to decide whether the matrix will be stored column-wise or row-wise. Here there is a huge difference with a sparse matrix where one can find many different storage schemes, see for example Section 2.1 *Storage Formats* in the SPARSKIT paper

http://www-users.cs.umn.edu/~saad/software/SPARSKIT/index.html

Each storage scheme has its own advantages and disadvantages. When we use some sparse solver, then the storage scheme is usually fixed and we must supply a matrix in some predefined format. However, such a format may not be the best to prepare the matrix or it well might be that the sparse matrix format in the code is already defined based on some other considerations. Hence in practice it is necessary anyway to convert a sparse matrix from one storage scheme to another.

I describe here a simple sparse matrix class that is well suited when the number of nonzero entries is not known in advance and memory should be dynamically allocated, for example assembly in finite elements of a global matrix from element matrices. In my code I will use it to read the matrix in the Matrix Market format and then the matrix will be transformed to the storage required by a sparse solver.

A sparse matrix could be considred as a list of triples: a nonzero entry, a row index, a column index. For example, one could imagine

struct entry{ double val; size_t i; size_t j;};

vector<struc> A:

It is easy to allocate entry memory but the problem here is a search for a particular element, for example `A(i, j)`

is not a trivial operation. Also it is necessary to pay attention that theoretically it is possible that such a structure will contain two entries with the same *i* and *j*.

After some trials, I have decided to employ

`vector<map<int, double> > A;`

where the matrix is supposed to be stored in the vector row-wise. That is, `A[i]`

returns a nonzero entries in the *i*-th row. `map<int, double>`

on the other hand is an efficient means to retrive the entry in the *i*-th row that corresponds to the *j*-th columns.

This means that `A[i][j]`

retrives the `A_ij`

and if there was no such entry it will be created. For example

`A[i][j] += val`

will always work correctly and efficiently, independently whether such an entry was already in the sparse matrix or not. When the matrix is ready it is possible to iterate through all nonzero entries as follows

for (int i = 0; i < A.size(); ++i) for (map<int, double>::const_iterator j = A[i].begin(); j != A[i].end(); ++j) // at this point: // i is the row index // (*j).first is the column index // (*j).second is value

The drawback here is that is requires some memory to store an array of trees. However, the matrix later will be converted to some more memory efficient scheme accepted by a sparse solver and the memory required to prepare the matrix will be freed before the linear solve.

The class `SparseMatrix`

is derived from `vector<map<int, double> >`

, so it iherits all operations mentioned above directly.

`class SparseMatrix : public vector< map<int, double> >`

{

Additionally the class has three members

`size_t nnz;`

size_t m;

bool IsSymmetric;

to keep the number of columns, the number of nonzero entries, and whether the matrix is stored completely (`IsSymmetric == false`

) or in the case of a symmetric matrix only upper part is available (`IsSymmetric == true`

). The functions of the class

`void makeNotSymmetric();`

size_t NRows() const {return size();}

size_t NCols() const {return m;}

void swap(SparseMatrix &y);

void clearMemory();

void read(istream &in);

void write(ostream &out);

are similar to those described for class `Matrix`

. The function `makeNotSymmetric`

converts the matrix from symmetric storage to the full storage.

The declaration of the class `SparseMatrix`

is in matrices.h and the definitions are in matrices.cpp.

## Previous

## Next

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

## Discussion

SparseMatrix : deallocating some elements

https://groups.google.com/d/topic/matrixprogramming/s0H8D6nY2VI/discussion

RSS