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

Introduction
Class Matrix

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