improve slicing lazily transposed sparse matrices: WIP

[jlapeyre]

This PR uses the code for getting slices in sparsevector.jl to get slices of sparse matrices
that have lazy adjoint and transpose wrappers. Before, the fallback code in abstractarray.jl created the
slices inefficiently.

Here is an example of a sparse matrix created with sprand(10^4 10^4, 0.01), which shows an increase in speed of 10 to 1000 x.
The test code.
Benchmark before the PR.
Benchmark after the PR

The changes should have no effect on the efficiency of slicing unwrapped, ie SparseMatrixCSC, matrices.

The new slicing calls adjoint and transpose recursively. I don't know of any cases where that is useful, although there may be. Consider for example a sparse matrix with elements that are dense matrices. If you want to interpret this as a block matrix, then transpose(S)[1, :] doesn't give you the correct result. (Neither does S[1,:].) Maybe restricting the recursion to structures that are intended to be used for linear algebra (eg BlockArray), rather than any ad-hoc array makes sense. The reason for this restriction would be to reduce complexity. I'm sure this has been discussed elsewhere, but I have not yet read about it.

There could be further optimizations. EDIT: The following has been done and commited to this branch.
For instance, check if the elements of the parent matrix are numbers (via dispatch), in which case it is not necessary to call transpose.

EDIT: Tests have been added to this branch This PR needs tests, there are none. I would not be surprised if some or all slicing of wrapped arrays is not covered by existing tests.

[jlapeyre]

Ref to discourse: https://discourse.julialang.org/t/strange-performance-with-adjoint-structures/13397

[jlapeyre]

Testing with code coverage showed that before this PR there were no tests of indexing sparse matrices wrapped with Transpose or Adjoint.

This commit adds tests for all of the new indexing methods.

[jlapeyre]

Commit #e5d2acd decreases the benchmark time consistently by about 10% for some test matrices. I think the compiler has enough information to do the same optimization without this commit. But, that optimization might not come soon.

[chriscoey]

nice! I'm looking forward to this being merged. is it still WIP?

[jlapeyre]

It is finished. I just needs to be reviewed. I changed existing indexing routines. But, I checked that the changes are optimized away, and perform exactly as before. Still someone should trigger nanosoldier.

[andreasnoack]

@nanosoldier runbenchmarks("sparse", vs = ":master")

[jlapeyre]

Where can I find the output from nanosoldier ? ... or maybe it failed to run.

[andreasnoack]

I think it was on strike when I tried last time but should be back on duty now so let's try again

@nanosoldier runbenchmarks("sparse", vs = ":master")

[andreasnoack]

Hm. It looks like I've been misinformed. It still isn't running.

[jlapeyre]

Sorry about your wasted time mucking with this. I doubt it will show a performance hit, but the modification really should be tested properly. I prefer to wait a bit rather than run them locally.

[andreasnoack]

@nanosoldier runbenchmarks("sparse", vs = ":master")

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @ararslan

[jlapeyre]

Looks like there may be a regression. And, reviewing my own benchmarks in the gist I see a regression:

# before PR
@btime a[1, :];
 25.862 μs (15 allocations: 4.56 KiB)

# after PR
@btime a[1, :];
  31.257 μs (16 allocations: 4.59 KiB)

These should be unchanged. But, in fact there is an additional allocation. I thought I took care of this before the PR. It can definitely be fixed by copying code. But, I want to avoid that, if possible. I will look into it...

[ViralBShah]

Let's try to get this in.

[ViralBShah]

Bump.

[jlapeyre]

I have four PRs to Julia in the works. I'll prioritize this one. But, at the moment, my time budget is pretty tight. Thanks for the reminder.

[ViralBShah]

Just checking to see if it is worth revisiting now.

[DilumAluthge]

We have moved the SparseArrays stdlib to an external repository.

Please open this PR against that repository: https://github.com/JuliaLang/SparseArrays.jl

Thank you!

[jlapeyre]

Huh. I completely forgot about this. I'll look into it.