[Rcpp-devel] Sparse matrices with RcppArmadillo

Hi Romain & Doug,

Two points to note:

(1)
When poking around internal Armadillo pointers and arrays for sparse
matrices, please take into account that the col_ptrs array is slightly
longer than specified by the Compressed Sparse Column (CSC) format.

However, this should not cause any problems when interfacing with
external libraries. If you use the set_size() function before writing
to the values[], row_indices[] and col_ptrs[] arrays, everything
should be okay.

The col_ptrs array is slightly longer in order to allow robust
iterators.  Specifically, instead of the col_ptrs array having a size
of ncols+1 (as specified by CSC), it has a size of ncols+2.  The extra
element is always set to std::numeric_limits<uword>::max().  As such,
be careful not to overwrite this magic number.

For more details within the code, see the function
SpMat<eT>::init(uword in_rows, uword in_cols), on line 3635 in
"SpMat_meat.hpp".

(2)
If you need to allocate memory for the values[] and row_indices[]
arrays, you need to use the memory::acquire_chunked() function.  This
is an alternative memory allocation mechanism designed to reduce the
burden of inserting or removing a non-zero element.  The rest of the
sparse matrix memory management code relies on this.  See the code in
file "memory.hpp" for more details.


On Sun, Dec 9, 2012 at 5:55 AM, Romain Francois

Yes. Conrad left a back door open so that we can add things **inside** the
SpMat template, similar to what we have done in dense matrix classes.

Going from there I was looking at the copy constructor of a SpMat, which
eventually is a call to init:

/**
 * Copy from another matrix.
 */
template<typename eT>
inline
void
SpMat<eT>::init(const SpMat<eT>& x)
  {
  arma_extra_debug_sigprint();

  // Ensure we are not initializing to ourselves.
  if (this != &x)
    {
    init(x.n_rows, x.n_cols);

    // values and row_indices may not be null.
    if (values != NULL)
      {
      memory::release(values);
      memory::release(row_indices);
      }

    access::rw(values)      = memory::acquire_chunked<eT> (x.n_nonzero + 1);
    access::rw(row_indices) = memory::acquire_chunked<uword>(x.n_nonzero +
1);

    // Now copy over the elements.
    arrayops::copy(access::rwp(values),      x.values,      x.n_nonzero +
1);
    arrayops::copy(access::rwp(row_indices), x.row_indices, x.n_nonzero +
1);
    arrayops::copy(access::rwp(col_ptrs),    x.col_ptrs,    x.n_cols + 1);

    access::rw(n_nonzero) = x.n_nonzero;
    }
  }

So it all looks very doable.

Romain