Skip to content

Commit 2b878f0

Browse files
authored
LinearAlgebra: Correct zero element in _generic_matvecmul! for block 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] ```
1 parent 4a4d850 commit 2b878f0

File tree

2 files changed

+16
-2
lines changed

2 files changed

+16
-2
lines changed

stdlib/LinearAlgebra/src/matmul.jl

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -779,7 +779,8 @@ function _generic_matvecmul!(C::AbstractVector, tA, A::AbstractVecOrMat, B::Abst
779779
else
780780
for k = 1:mA
781781
aoffs = (k-1)*Astride
782-
s = zero(A[aoffs + 1]*B[1] + A[aoffs + 1]*B[1])
782+
firstterm = transpose(A[aoffs + 1])*B[1]
783+
s = zero(firstterm + firstterm)
783784
for i = 1:nA
784785
s += transpose(A[aoffs+i]) * B[i]
785786
end
@@ -794,7 +795,8 @@ function _generic_matvecmul!(C::AbstractVector, tA, A::AbstractVecOrMat, B::Abst
794795
else
795796
for k = 1:mA
796797
aoffs = (k-1)*Astride
797-
s = zero(A[aoffs + 1]*B[1] + A[aoffs + 1]*B[1])
798+
firstterm = A[aoffs + 1]'B[1]
799+
s = zero(firstterm + firstterm)
798800
for i = 1:nA
799801
s += A[aoffs + i]'B[i]
800802
end

stdlib/LinearAlgebra/test/matmul.jl

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -216,6 +216,18 @@ end
216216
end
217217
end
218218

219+
@testset "generic_matvecmul for vectors of matrices" begin
220+
x = [1 2 3; 4 5 6]
221+
A = reshape([x,2x,3x,4x],2,2)
222+
b = [x, 2x]
223+
for f in (adjoint, transpose)
224+
c = f(A) * b
225+
for i in eachindex(c)
226+
@test c[i] == sum(f(A)[i, j] * b[j] for j in eachindex(b))
227+
end
228+
end
229+
end
230+
219231
@testset "generic_matmatmul for matrices of vectors" begin
220232
B = Matrix{Vector{Int}}(undef, 2, 2)
221233
B[1, 1] = [1, 2]

0 commit comments

Comments
 (0)