Skip to content

Commit

Permalink
LinearAlgebra: Correct zero element in _generic_matvecmul! for bloc…
Browse files Browse the repository at this point in the history
…k adj/trans (#54151)

Fixes the following issue on master, where the zero element is computed
incorrectly (but subsequent terms are computed correctly):
```julia
julia> using LinearAlgebra

julia> x = [1 2 3; 4 5 6];

julia> A = reshape([x,2x,3x,4x],2,2);

julia> b = fill(x, 2);

julia> A' * b
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(2),Base.OneTo(3)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))
Stacktrace:
  [1] _generic_matmatmul!(C::Matrix{Int64}, A::Matrix{Int64}, B::Matrix{Int64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
    @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:849
  [2] generic_matmatmul!
    @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:834 [inlined]
  [3] _mul!
    @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:287 [inlined]
  [4] mul!
    @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:285 [inlined]
  [5] mul!
    @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:253 [inlined]
  [6] *
    @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:114 [inlined]
  [7] _generic_matvecmul!(C::Vector{Matrix{…}}, tA::Char, A::Matrix{Matrix{…}}, B::Vector{Matrix{…}}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
    @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:797
  [8] generic_matvecmul!
    @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:755 [inlined]
  [9] _mul!
    @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:73 [inlined]
 [10] mul!
    @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:70 [inlined]
 [11] mul!
    @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:253 [inlined]
 [12] *(A::Adjoint{Adjoint{Int64, Matrix{Int64}}, Matrix{Matrix{Int64}}}, x::Vector{Matrix{Int64}})
    @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:60
 [13] top-level scope
    @ REPL[10]:1
Some type information was truncated. Use `show(err)` to see complete types.
```
After this PR,
```julia
julia> A' * b
2-element Vector{Matrix{Int64}}:
 [51 66 81; 66 87 108; 81 108 135]
 [119 154 189; 154 203 252; 189 252 315]
```

(cherry picked from commit 2b878f0)
  • Loading branch information
jishnub authored and KristofferC committed May 8, 2024
1 parent 4d21550 commit 280eac4
Show file tree
Hide file tree
Showing 2 changed files with 16 additions and 2 deletions.
6 changes: 4 additions & 2 deletions stdlib/LinearAlgebra/src/matmul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -711,7 +711,8 @@ function _generic_matvecmul!(C::AbstractVector, tA, A::AbstractVecOrMat, B::Abst
else
for k = 1:mA
aoffs = (k-1)*Astride
s = zero(A[aoffs + 1]*B[1] + A[aoffs + 1]*B[1])
firstterm = transpose(A[aoffs + 1])*B[1]
s = zero(firstterm + firstterm)
for i = 1:nA
s += transpose(A[aoffs+i]) * B[i]
end
Expand All @@ -726,7 +727,8 @@ function _generic_matvecmul!(C::AbstractVector, tA, A::AbstractVecOrMat, B::Abst
else
for k = 1:mA
aoffs = (k-1)*Astride
s = zero(A[aoffs + 1]*B[1] + A[aoffs + 1]*B[1])
firstterm = A[aoffs + 1]'B[1]
s = zero(firstterm + firstterm)
for i = 1:nA
s += A[aoffs + i]'B[i]
end
Expand Down
12 changes: 12 additions & 0 deletions stdlib/LinearAlgebra/test/matmul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,18 @@ end
end
end

@testset "generic_matvecmul for vectors of matrices" begin
x = [1 2 3; 4 5 6]
A = reshape([x,2x,3x,4x],2,2)
b = [x, 2x]
for f in (adjoint, transpose)
c = f(A) * b
for i in eachindex(c)
@test c[i] == sum(f(A)[i, j] * b[j] for j in eachindex(b))
end
end
end

@testset "generic_matmatmul for matrices of vectors" begin
B = Matrix{Vector{Int}}(undef, 2, 2)
B[1, 1] = [1, 2]
Expand Down

0 comments on commit 280eac4

Please sign in to comment.