[Octave-patch-tracker] [patch #8693] Reimplement operator * (const SparseMatrix& m, const SparseMatrix& a) to use SuiteSparse

[anonymous]

URL:
  <http://savannah.gnu.org/patch/?8693>

                 Summary: Reimplement operator * (const SparseMatrix& m, const
SparseMatrix& a) to use SuiteSparse
                 Project: GNU Octave
            Submitted by: None
            Submitted on: Tue 23 Jun 2015 08:30:45 PM UTC
                Category: None
                Priority: 5 - Normal
                  Status: None
                 Privacy: Public
             Assigned to: None
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any

    _______________________________________________________

Details:

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). 



    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Tue 23 Jun 2015 08:30:45 PM UTC  Name: suitesparse_mul.patch  Size: 7kB 
 By: None

<http://savannah.gnu.org/patch/download.php?file_id=34300>

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?8693>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/