Skip to content

[BUG] Block sparse matvec is broken #199

@mtfishman

Description

@mtfishman

For example:

julia> using BlockSparseArrays

julia> r = [2, 2];

julia> a1 = BlockSparseArray{Float64}(undef, r, r);

julia> a2 = BlockSparseArray{Float64}(undef, r);

julia> a1 * a2
#...

julia> err
1-element ExceptionStack:
Not implemented.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:44
  [2] materialize!(m::ArrayLayouts.MulAdd{BlockArrays.BlockLayout{SparseArraysBase.SparseLayout, ArrayLayouts.DenseColumnMajor}, BlockArrays.BlockLayout{SparseArraysBase.SparseLayout, ArrayLayouts.DenseColumnMajor}, BlockArrays.BlockLayout{SparseArraysBase.SparseLayout, ArrayLayouts.DenseColumnMajor}, Float64, BlockSparseMatrix{Float64, Matrix{Float64}, SparseArraysBase.SparseMatrixDOK{Matrix{Float64}, BlockSparseArrays.ZeroBlocks{2, Matrix{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}, BlockSparseVector{Float64, Vector{Float64}, SparseArraysBase.SparseVectorDOK{Vector{Float64}, BlockSparseArrays.ZeroBlocks{1, Vector{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}, BlockSparseVector{Float64, Vector{Float64}, SparseArraysBase.SparseVectorDOK{Vector{Float64}, BlockSparseArrays.ZeroBlocks{1, Vector{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}})
    @ BlockSparseArrays ~/.julia/packages/BlockSparseArrays/gEM83/src/blocksparsearrayinterface/arraylayouts.jl:31
  [3] muladd!(α::Float64, A::BlockSparseMatrix{Float64, Matrix{Float64}, SparseArraysBase.SparseMatrixDOK{Matrix{Float64}, BlockSparseArrays.ZeroBlocks{2, Matrix{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}, B::BlockSparseVector{Float64, Vector{Float64}, SparseArraysBase.SparseVectorDOK{Vector{Float64}, BlockSparseArrays.ZeroBlocks{1, Vector{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}, β::Float64, C::BlockSparseVector{Float64, Vector{Float64}, SparseArraysBase.SparseVectorDOK{Vector{Float64}, BlockSparseArrays.ZeroBlocks{1, Vector{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}; kw::@Kwargs{Czero::Bool})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/RCHQ8/src/muladd.jl:75
  [4] copyto!
    @ ~/.julia/packages/ArrayLayouts/RCHQ8/src/muladd.jl:82 [inlined]
  [5] copy(M::ArrayLayouts.MulAdd{BlockArrays.BlockLayout{SparseArraysBase.SparseLayout, ArrayLayouts.DenseColumnMajor}, BlockArrays.BlockLayout{SparseArraysBase.SparseLayout, ArrayLayouts.DenseColumnMajor}, ArrayLayouts.ZerosLayout, Float64, BlockSparseMatrix{Float64, Matrix{Float64}, SparseArraysBase.SparseMatrixDOK{Matrix{Float64}, BlockSparseArrays.ZeroBlocks{2, Matrix{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}, BlockSparseVector{Float64, Vector{Float64}, SparseArraysBase.SparseVectorDOK{Vector{Float64}, BlockSparseArrays.ZeroBlocks{1, Vector{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}, FillArrays.Zeros{Float64, 1, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/RCHQ8/src/muladd.jl:77
  [6] copy
    @ ~/.julia/packages/ArrayLayouts/RCHQ8/src/mul.jl:140 [inlined]
  [7] materialize
    @ ~/.julia/packages/ArrayLayouts/RCHQ8/src/mul.jl:137 [inlined]
  [8] mul
    @ ~/.julia/packages/ArrayLayouts/RCHQ8/src/mul.jl:138 [inlined]
  [9] *(A::BlockSparseMatrix{Float64, Matrix{Float64}, SparseArraysBase.SparseMatrixDOK{Matrix{Float64}, BlockSparseArrays.ZeroBlocks{2, Matrix{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}, B::BlockSparseVector{Float64, Vector{Float64}, SparseArraysBase.SparseVectorDOK{Vector{Float64}, BlockSparseArrays.ZeroBlocks{1, Vector{Float64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/RCHQ8/src/mul.jl:226
 [10] top-level scope
    @ REPL[49]:1
 [11] top-level scope
    @ REPL:1

A workaround is adding a trivial dimension:

julia> a2 = BlockSparseArray{Float64}(undef, r, [1]);

julia> a1 * a2
2×1-blocked 4×1 BlockSparseMatrix{Float64, Matrix{Float64}, , }:
   
   
 ───
   
   

See #197.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions