patchGNU Octave - Patches: patch #8693, Reimplement operator * (const...


patch #8693: Reimplement operator * (const SparseMatrix& m, const SparseMatrix& a) to use SuiteSparse

Submitted by:  None
Submitted on:  Tue 23 Jun 2015 08:30:45 PM UTC  
Category: Core : otherPriority: 5 - Normal
Status: In ProgressPrivacy: Public
Assigned to: NoneOriginator Email: -unavailable-
Open/Closed: Open

Add a New Comment(Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


Sun 10 Jul 2016 12:29:21 AM UTC, comment #4:

The benefit of this patch is a performance improvement on sparse multiplication. I think it's worth considering if someone who is familiar with the routines involved, maybe has some sparse code they can run against this, can take a look and test. I don't use sparse matrices and I'm not familiar with the libraries at all.

Mike Miller <mtmiller>
Project Administrator
Wed 06 Jul 2016 01:58:23 PM UTC, comment #3:

Mike, do you think this should be considered for 4.2.0? If so, do you know who could comment on the design choices, which you declined to comment on?

Lachlan Andrew <lachlan>
Project Member
Sun 13 Dec 2015 08:04:42 PM UTC, comment #2:

I've reviewed this patch for coding style, commit message guidelines, and verified that it applies cleanly, compiles, and make check runs with no new failures.

I have not reviewed the design choices or the basic functionality or performance changes that this patch addresses.

Revised patch attached with the following changes:

  • author's name in the commit header
  • commit message edited to be consistent with our practice
  • whitespace cleanup
  • TODO changed to FIXME
  • the last line of the function changed to compile for me (and reason for the const_cast unnecessary I think)

There's not much more reviewing I can bring to this patch, so not assigning to myself. Someone with more familiarity with the design of the sparse array classes and with the SuiteSparse library will have to validate the rest.

(file #35730)

Mike Miller <mtmiller>
Project Administrator
Fri 24 Jul 2015 12:30:05 PM UTC, comment #1:

As suggested on the #octave IRC channel, I rewrote the submission as a patch, and added comments through the code.

(file #34502)

Daniele Calandriello <danielec>
Tue 23 Jun 2015 08:30:45 PM UTC, original submission:


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



(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):

Attached Files
file #35730:  revised-8693.patch added by mtmiller (10KiB - text/x-diff)
file #34502:  patch_suitesparse_mul#8693.cset added by danielec (10KiB - application/octet-stream)
file #34300:  suitesparse_mul.patch added by None (7KiB - text/x-patch)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -unavailable- added by lachlan (Posted a comment)
  • -unavailable- added by mtmiller (Updated the item)
  • -unavailable- added by danielec (Updated the item)

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only logged-in users can vote.


    Please enter the title of George Orwell's famous dystopian book (it's a date):



    Follow 5 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sun 10 Jul 2016 12:29:21 AM UTCmtmillerCategoryNone=>Core : other
    Sun 13 Dec 2015 08:04:42 PM UTCmtmillerAttached File-=>Added revised-8693.patch, #35730
      StatusNone=>In Progress
    Fri 24 Jul 2015 12:30:05 PM UTCdanielecAttached File-=>Added patch_suitesparse_mul#8693.cset, #34502
    Tue 23 Jun 2015 08:30:45 PM UTCNoneAttached File-=>Added suitesparse_mul.patch, #34300

    Back to the top

    Powered by Savane 3.1-cleanup1