BlockArray to Sparse quite slow

[rveltz]

Hi,

I noticed the following is quite slow

using LinearAlgebra, SparseArrays, BlockArrays

M = 30
N = 500
J = BlockArray(spzeros(M * N, M * N), N * ones(Int64,M),  N * ones(Int64,M))

for ii=2:M
    setblock!(J, sprand(N,N,1/N), ii, ii)
    setblock!(J, sprand(N,N,1/N), ii,ii-1)
end

J1 = @time sparse(J) #5seconds

whereas

function createSP(J::BlockArray)
    nl, nc = size(J.blocks)
    N = size(J[Block(1,1)])[1]
    res = spzeros(N,N*nc)
    for i=1:nl
        line = J[Block(i,1)]
        for j=2:nc
            line = hcat(line,J[Block(i,j)])
        end
        res = vcat(res,line)
    end
    return res[N+1:end,:]
end

Jb = @time createSP(J)

takes 0.23s.

1. Am I missing a point?
2. Any chance to put createSP in your package?

Thank you,

Best regards

[dlfivefifty]

sparse(::BlockArray) is just calling the default, which spends a lot of time in findall(x -> x != 0, M). It's possible a faster findall would accomplish the same feat, otherwise, can you file a PR with your createSP implemented via overloading SparseMatrixCSC{Tv,Ti}(M::AbstractBlockArray) where {Tv,Ti} which sparse eventually calls?

There's a question whether sparse(::AbstractBlockArray) should return a block array itself with sparse storage, so as to not lose the block information?

[rveltz]

We should limit the type to BlockArray{Tv,2,SparseMatrixCSC{Tv,Ti}}, don't you agree?

I am not sure I am so skilled as to do a pull request ;)

I wanted to check if I missed a point. Also, my function is geared toward blocks with same size.

[dlfivefifty]

We should limit the type to BlockArray{Tv,2,SparseMatrixCSC{Tv,Ti}}, don't you agree?

No, as other types have nice sparse representations. Calling sparse(J[Block(i,j)]) should tackle the general case.

But it may be better to create the i, j, z directly, say by modifying sparse(::BandedBlockBandedMatrix): https://github.com/JuliaMatrices/BlockBandedMatrices.jl/blob/e19c22fc9c9e5399976c24aa06a6871bf0133841/src/BandedBlockBandedMatrix.jl#L254

I am not sure I am so skilled as to do a pull request ;)

In open-source software, people maintaining packages have limited time and so are unlikely to add new features that they don't directly benefit from. So if you really want this you'll probably have to do the PR yourself. Given the current state of Julia packages I'd consider being able to make a PR an essential skill to have.

[rveltz]

Which file should I put the function in?

[dlfivefifty]

If it only works with BlockArray, then probably src/blockarray.jl, otherwise src.abstractblockarray.jl

[briochemc]

Did a PR eventually come out of this?

[dlfivefifty]

I don't believe so, if you are keen to do one I'm sure it will be appreciated.

[briochemc]

I was in need of the same thing and thought: What about

Jcat = @time reduce(hcat, reduce(vcat, J[Block(i, j)] for i in 1:blocksize(J, 1)) for j in 1:blocksize(J, 2))

No idea if this is not good in some cases but I got similar times as @rveltz

julia> Jb = @time createSP(J)
  0.139867 seconds (95.83 k allocations: 132.349 MiB, 51.51% compilation time)
15000×15000 SparseMatrixCSC{Float64, Int64} with 28666 stored entries:
⎡⣤⣤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠻⢿⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠈⠛⣿⣶⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠘⠻⢿⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠈⠛⣿⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠈⠹⣿⣧⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠛⣿⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠹⣿⣧⣀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⢿⣦⡄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣿⣧⣀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⢿⣦⡄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣿⣶⎦

julia> Jcat = @time reduce(hcat, reduce(vcat, J[Block(i, j)] for i in 1:blocksize(J, 1)) for j in 1:blocksize(J, 2))
  0.141077 seconds (282.44 k allocations: 37.092 MiB, 87.99% compilation time)
15000×15000 SparseMatrixCSC{Float64, Int64} with 28666 stored entries:
⎡⣤⣤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠻⢿⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠈⠛⣿⣶⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠘⠻⢿⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠈⠛⣿⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠈⠹⣿⣧⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠛⣿⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠹⣿⣧⣀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⢿⣦⡄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣿⣧⣀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⢿⣦⡄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣿⣶⎦

so maybe this is good enough?