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


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

Submitter:  None
Submitted:  Tue 23 Jun 2015 08:30:45 PM UTC
Category:  Core : new feature Priority:  5 - Normal
Status:  None Privacy:  Public
Assigned to:  None Originator Email:  -email is unavailable-
Open/Closed:  Open
* Mandatory Fields

Add a New Comment Rich Markup

Sat 18 Jan 2020 11:05:34 AM UTC, comment #5: 

Review for Octave 6.

No progress on this item for 4 years.  Still interesting.

Kai Torben Ohlhus <siko1056>
Group Member
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>
Group 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>
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>
Group 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 Files:

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
  • -email is unavailable- added by siko1056 (Posted a comment)
  • -email is unavailable- added by lachlan (Posted a comment)
  • -email is unavailable- added by danielec (Updated the item)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.


    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-01-25 mtmiller Carbon-CopyRemoved 80942 -
    2020-01-18 siko1056 CategoryCore : other Core : new feature
        StatusIn Progress None
    2016-07-10 mtmiller CategoryNone Core : other
    2015-12-13 mtmiller Attached File- Added revised-8693.patch, #35730
        StatusNone In Progress
    2015-07-24 danielec Attached File- Added patch_suitesparse_mul#8693.cset, #34502
    2015-06-23 None Attached File- Added suitesparse_mul.patch, #34300

    Back to the top

    Powered by Savane 3.12