mul! with two BandedMatrices that share data through a view throws an ArgumentError

[MikaelSlevinsky]

mul! can work with nonintersecting views of dense matrices:

julia> using LinearAlgebra

julia> A = zeros(10, 10); A[1:5, 1:5] .= randn(5, 5);

julia> B = randn(5,5);

julia> A1 = view(A, 1:5, 1:5);

julia> A2 = view(A, 1:5, 6:10);

julia> using LinearAlgebra

julia> A = zeros(5, 10); A[1:5, 1:5] .= randn(5, 5);

julia> B = randn(5,5);

julia> A1 = view(A, 1:5, 1:5);

julia> A2 = view(A, 1:5, 6:10);

julia> mul!(A2, B, A1)
5×5 view(::Matrix{Float64}, 1:5, 6:10) with eltype Float64:
 -2.65047   1.63847   -0.972346  -0.803627    2.03012
 -0.817101  1.09273   -1.04688    0.41023    -1.15711
 -3.21165   0.721003  -5.99989   -2.00797     7.1199
 -1.67402   1.31855   -4.07428   -0.0920713   1.73292
 -1.01159   2.77983    1.81222    0.589318   -2.18872

Success!

With banded matrices, this same kind of operation doesn't work:

julia> using BandedMatrices

julia> ADATA = randn(3, 10);

julia> A1 = BandedMatrices._BandedMatrix(view(ADATA, 1:3, 1:5), 5, 1, 1);

julia> A2 = BandedMatrices._BandedMatrix(view(ADATA, 1:3, 6:10), 5, 1, 1);

julia> mul!(A2, B, A1)
ERROR: ArgumentError: an array of type `BandedMatrix` shares memory with another argument
and must make a preventative copy of itself in order to maintain consistent semantics,
but `copy(::BandedMatrix{Float64, SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, Base.OneTo{Int64}})` returns a new array of type `BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}`.
To fix, implement:
    `Base.unaliascopy(A::BandedMatrix)::typeof(A)`
Stacktrace:
 [1] _unaliascopy(A::BandedMatrix{Float64, SubArray{…}, Base.OneTo{…}}, C::BandedMatrix{Float64, Matrix{…}, Base.OneTo{…}})
   @ Base ./abstractarray.jl:1500
 [2] unaliascopy(A::BandedMatrix{Float64, SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{…}, UnitRange{…}}, false}, Base.OneTo{Int64}})
   @ Base ./abstractarray.jl:1498
 [3] unalias
   @ ./abstractarray.jl:1481 [inlined]
 [4] copyto!
   @ ~/.julia/packages/ArrayLayouts/B2wRU/src/muladd.jl:82 [inlined]
 [5] copyto!
   @ ~/.julia/packages/ArrayLayouts/B2wRU/src/mul.jl:142 [inlined]
 [6] mul!(dest::BandedMatrix{Float64, SubArray{…}, Base.OneTo{…}}, A::Matrix{Float64}, B::BandedMatrix{Float64, SubArray{…}, Base.OneTo{…}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/B2wRU/src/mul.jl:143
 [7] mul!(dest::BandedMatrix{Float64, SubArray{…}, Base.OneTo{…}}, A::Matrix{Float64}, B::BandedMatrix{Float64, SubArray{…}, Base.OneTo{…}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/B2wRU/src/mul.jl:212
 [8] top-level scope
   @ REPL[31]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> 

[MikaelSlevinsky]

Of course, B*A1 above would not be banded, but the error is the same when B is a banded matrix with (0,0) bandwidths:

julia> B = BandedMatrices._BandedMatrix(randn(1, 5), 5, 0, 0);

julia> mul!(A2, B, A1)
ERROR: ArgumentError: an array of type `BandedMatrix` shares memory with another argument
and must make a preventative copy of itself in order to maintain consistent semantics,
but `copy(::BandedMatrix{Float64, SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, Base.OneTo{Int64}})` returns a new array of type `BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}`.
To fix, implement:
    `Base.unaliascopy(A::BandedMatrix)::typeof(A)`
Stacktrace:
 [1] _unaliascopy(A::BandedMatrix{Float64, SubArray{…}, Base.OneTo{…}}, C::BandedMatrix{Float64, Matrix{…}, Base.OneTo{…}})
   @ Base ./abstractarray.jl:1500
 [2] unaliascopy(A::BandedMatrix{Float64, SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{…}, UnitRange{…}}, false}, Base.OneTo{Int64}})
   @ Base ./abstractarray.jl:1498
 [3] unalias
   @ ./abstractarray.jl:1481 [inlined]
 [4] copyto!
   @ ~/.julia/packages/ArrayLayouts/B2wRU/src/muladd.jl:82 [inlined]
 [5] copyto!
   @ ~/.julia/packages/ArrayLayouts/B2wRU/src/mul.jl:142 [inlined]
 [6] mul!(dest::BandedMatrix{…}, A::BandedMatrix{…}, B::BandedMatrix{…})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/B2wRU/src/mul.jl:143
 [7] mul!(dest::BandedMatrix{…}, A::BandedMatrix{…}, B::BandedMatrix{…})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/B2wRU/src/mul.jl:214
 [8] top-level scope
   @ REPL[38]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> 

and, interestingly, if B = Diagonal(randn(5)); then there's a stack overflow 🤯.

[jishnub]

Which version are you using? I obtain an error in the first case, but a different one:

julia> using BandedMatrices, LinearAlgebra

julia> B = randn(5,5);

julia> ADATA = randn(3, 10);

julia> A1 = BandedMatrices._BandedMatrix(view(ADATA, 1:3, 1:5), 5, 1, 1);

julia> A2 = BandedMatrices._BandedMatrix(view(ADATA, 1:3, 6:10), 5, 1, 1);

julia> mul!(A2, B, A1)
ERROR: BandError: attempt to access SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false} with bandwidths (2, 4) at band -2
Stacktrace:
  [1] banded_setindex!
    @ ~/.julia/packages/BandedMatrices/8wPaL/src/banded/BandedMatrix.jl:613 [inlined]
  [2] setindex!
    @ ~/.julia/packages/BandedMatrices/8wPaL/src/banded/BandedMatrix.jl:622 [inlined]
  [3] _default_blasmul_loop!
    @ ~/.julia/packages/ArrayLayouts/B2wRU/src/muladd.jl:169 [inlined]
  [4] default_blasmul!(α::Float64, A::Matrix{…}, B::BandedMatrix{…}, β::Bool, C::BandedMatrix{…})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/B2wRU/src/muladd.jl:188
  [5] materialize!(M::ArrayLayouts.MulAdd{…})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/B2wRU/src/muladd.jl:264
  [6] muladd!
    @ ~/.julia/packages/ArrayLayouts/B2wRU/src/muladd.jl:75 [inlined]
  [7] copyto!
    @ ~/.julia/packages/ArrayLayouts/B2wRU/src/muladd.jl:82 [inlined]
  [8] copyto!
    @ ~/.julia/packages/ArrayLayouts/B2wRU/src/mul.jl:142 [inlined]
  [9] mul!(dest::BandedMatrix{…}, A::Matrix{…}, B::BandedMatrix{…})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/B2wRU/src/mul.jl:143
 [10] mul!(dest::BandedMatrix{…}, A::Matrix{…}, B::BandedMatrix{…})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/B2wRU/src/mul.jl:212
 [11] top-level scope
    @ REPL[20]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> B = BandedMatrices._BandedMatrix(randn(1, 5), 5, 0, 0);

julia> mul!(A2, B, A1)
5×5 BandedMatrix{Float64} with bandwidths (1, 1) with data 3×5 view(::Matrix{Float64}, 1:3, 6:10) with eltype Float64:
 0.250817  -0.0113354     ⋅           ⋅         ⋅ 
 0.17432    0.00478642   0.0268603    ⋅         ⋅ 
  ⋅        -0.880198    -0.0148332  -0.208889   ⋅ 
  ⋅          ⋅          -0.293474   -0.691504  1.9991
  ⋅          ⋅            ⋅          0.269405  0.607511

julia> B = Diagonal(rand(5));

julia> mul!(A2, B, A1)
5×5 BandedMatrix{Float64} with bandwidths (1, 1) with data 3×5 view(::Matrix{Float64}, 1:3, 6:10) with eltype Float64:
 -0.113272   0.364919    ⋅           ⋅         ⋅ 
  0.119846  -1.87337   -0.0520645    ⋅         ⋅ 
   ⋅         0.148692   3.24714    -0.31883    ⋅ 
   ⋅          ⋅        -0.240135    0.106178  0.153283
   ⋅          ⋅          ⋅         -0.875653  0.930174

julia> VERSION
v"1.11.5"

(jl_mqKB7j) pkg> st BandedMatrices
Status `/tmp/jl_mqKB7j/Project.toml`
  [aae01518] BandedMatrices v1.9.3