Hello,

I was browsing my profiling data for an oct file using SparseMatrix multiplication, and I found out that octave's sparse-sparse matrix multiplication is using an original implementation.

To make a small comparison, I rewrote the thing to use SuiteSparse _multiply. My naive benchmark show a performance improvement around 20%, but they are very artificial and I would not take them too seriously. On the other hand, as jwe mentioned on the IRC channel, using a library implementation can change the maintenance efforts going forward, so maybe this patch might be useful even if it does not improve performance.

As for the technical aspects of the patch itself, SPARSE_SPARSE_MUL is currently implemented as a macro.

It is used for double-double double-complex complex-double and complex-complex sparse matrix multiplication, my patch only replaces double-double implementation, but it can be extended without too much effort I believe (more on this later).

I copied without modifications all the initial checks and the logic for all corner cases. The reimplementation only covers the last case when an _expensive_ multiplication takes place. It simply passes the existing data to SuiteSparse, relying on the the fact that the arguments to _multiply are const, and no modification will take place.

The main problem is processing the result. I first build a Sparse<double> object to directly access the inner sparse_rep. I _transfer_ownership_ (I cannot think of a better naming) of the pointers returned by SuiteSparse to avoid copies. SuiteSparse uses compressed column space, but (unlike octave and matlab) does not keep the rows ordered. I need to reorder them, but also reorder the weights.

<small detour>
Sparse<double> exposes a method sort that seems to do exactly this, but simply calling it results in code that fails around 20 tests. Manual checking of tests in ichol such as

A2 = gallery ("poisson", 30);
opts.type = "nofill";
opts.michol = "off";
L = ichol (A2, opts);
assert (norm (A2 - L*L', "fro") / norm (A2, "fro"), 0.0893, 1e-4)

shows that there is some problem with the reordering. I am not submitting a bug report because it might as well be my fault for improperly calling sort on an unstable object.
</small detour>

Since I cannot use Sparse<double>::sort(), I sort the rows using octave_sort<octave_idx_type> but this solution is suboptimal. In particular, I need to sort the rows but also permute the weights, and I can only permute indices and then permute the weights. A further templatization of octave_sort<octave_idx_type>::sort (T *data, octave_idx_type *idx, octave_idx_type nel, Comp comp) might help with this. Other options are (on my side) to understand better Sparse<double>::sort() and the phantomatic sparse-sort from octave/liboctave/util/sparse-sort.h
Another small performance improvement could be achieved by not passing the size of the matrix to the Sparse<double> constructor (this allocates an n-long vector). It is very minor, and it fails a bunch of tests. If anyone is interested I can explain this a little bit more and add some debugging results.

Apart from this, the implementation works, and passes all tests.
 
As for extending this to complex multiplication, there are a few step to take. On a surface level it looks like you could simply template my implementation (with many more corner cases), and then call cs_ci_multiply or cs_cl_multiply instead of cs_di_multiply and cs_dl_multiply. I am not an expert in SuiteSparse (or complex analysis) and I might be missing something. For example one can look at SuiteSparse sources for inspiration. In suitesparse_4.2.1.orig.tar.gz the folder SuiteSparse/MATLAB_Tools/SSMULT contains sources for sparse-sparse multiplication mex replacements. These can be adapted, but the effort is non-trivial (I think) to transition them to octave types (e.g. mwIndex to octave_idx_type).