A.1.6.1 Array and Sparse Class Differences

The number of elements in a sparse matrix is considered to be the number of nonzero elements, rather than the product of the dimensions. Therefore,

SparseMatrix sm;
…
int nnz = sm.nelem ();

returns the number of nonzero elements (like the interpreter function nnz). If the user really requires the number of elements in the matrix, including the nonzero elements, they should use numel rather than nelem. Note that for very large matrices, where the product of the two dimensions is larger than the representation of an unsigned int, numel can overflow. An example is speye (1e6) which will create a matrix with a million rows and columns, but only a million nonzero elements. In this case, the number of rows multiplied by the number of columns is more than two hundred times the maximum value that can be represented by an unsigned 32-bit int. The use of numel should, therefore, be avoided unless it is known that it will not overflow.

Extreme care is also required when using the elem method or the () operator which perform essentially the same function. The reason is that if a sparse object is non-const, then Octave will assume that a request for a zero element in a sparse matrix is in fact a request to create this element so it can be filled. Therefore, a piece of code like

SparseMatrix sm;
…
for (int j = 0; j < nc; j++)
  for (int i = 0; i < nr; i++)
    std::cerr << " (" << i << "," << j << "): " << sm(i,j) << "\n";

is a great way of turning a sparse matrix into a dense one, and a very slow way at that since it reallocates the sparse object for each zero element in the matrix.

A simple way of preventing the above from happening is to create a temporary constant version of the sparse matrix. Note that only the container for the sparse matrix will be copied, while the actual representation of the data will be shared between the two versions of the sparse matrix; This is not a costly operation. The example above, re-written to prevent sparse-to-dense conversion, is

SparseMatrix sm;
…
const SparseMatrix tmp (sm);
for (int j = 0; j < nc; j++)
  for (int i = 0; i < nr; i++)
    std::cerr << " (" << i << "," << j << "): " << tmp(i,j) << "\n";

Finally, because the sparse types aren’t represented by a contiguous block of memory, the fortran_vec method of Array<T> is not available. It is, however, replaced by three separate methods ridx, cidx, and data, that access the raw compressed column format that Octave sparse matrices are stored in. These methods can be used in a manner similar to elem to allow the matrix to be accessed or filled. However, it is up to the user to respect the sparse matrix compressed column format or the matrix will become corrupted.