Skip to content

Commit

Permalink
Fix generic_matvecmul! for cases where zero(eltype(C)) is not defined.
Browse files Browse the repository at this point in the history 
Fixes #11978.

Fix type instability in generic_matmatmul! when addition changes the
element type, e.g. for Bool.
  • Loading branch information
@andreasnoack
andreasnoack committed Aug 21, 2015
1 parent da7cc55 commit 351e8af
Show file tree
Hide file tree
Showing 2 changed files with 85 additions and 24 deletions.
98 changes: 74 additions & 24 deletions base/linalg/matmul.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -371,14 +371,17 @@ function generic_matvecmul!{T,S,R}(C::AbstractVector{R}, tA, A::AbstractVecOrMat
if mA != length(C)
throw(DimensionMismatch("result C has length $(length(C)), needs length $mA"))
end
z = zero(R)

Astride = size(A, 1)

if tA == 'T' # fastest case
for k = 1:mA
aoffs = (k-1)*Astride
s = z
if mB == 0
s = zero(R)
else
s = zero(A[aoffs + 1]*B[1] + A[aoffs + 1]*B[1])
end
for i = 1:nA
s += A[aoffs+i].'B[i]
end
Expand All @@ -387,19 +390,29 @@ function generic_matvecmul!{T,S,R}(C::AbstractVector{R}, tA, A::AbstractVecOrMat
elseif tA == 'C'
for k = 1:mA
aoffs = (k-1)*Astride
s = z
if mB == 0
s = zero(R)
else
s = zero(A[aoffs + 1]*B[1] + A[aoffs + 1]*B[1])
end
for i = 1:nA
s += A[aoffs+i]'B[i]
s += A[aoffs + i]'B[i]
end
C[k] = s
end
else # tA == 'N'
fill!(C, z)
for i = 1:mA
if mB == 0
C[i] = zero(R)
else
C[i] = zero(A[i]*B[1] + A[i]*B[1])
end
end
for k = 1:mB
aoffs = (k-1)*Astride
b = B[k]
for i = 1:mA
C[i] += A[aoffs+i] * b
C[i] += A[aoffs + i] * b
end
end
end
Expand Down Expand Up @@ -490,27 +503,40 @@ function generic_matmatmul!{T,S,R}(C::AbstractVecOrMat{R}, tA, tB, A::AbstractVe
end
else
# Multiplication for non-plain-data uses the naive algorithm

if tA == 'N'
if tB == 'N'
for i = 1:mA, j = 1:nB
Ctmp = A[i, 1]*B[1, j]
for k = 2:nA
if length(A) == 0 || length(B) == 0
Ctmp = zero(R)
else
Ctmp = zero(A[i, 1]*B[1, j] + A[i, 1]*B[1, j])
end
for k = 1:nA
Ctmp += A[i, k]*B[k, j]
end
C[i,j] = Ctmp
end
elseif tB == 'T'
for i = 1:mA, j = 1:nB
Ctmp = A[i, 1]*B[j, 1].'
for k = 2:nA
if length(A) == 0 || length(B) == 0
Ctmp = zero(R)
else
Ctmp = zero(A[i, 1]*B[j, 1] + A[i, 1]*B[j, 1])
end
for k = 1:nA
Ctmp += A[i, k]*B[j, k].'
end
C[i,j] = Ctmp
end
else
for i = 1:mA, j = 1:nB
Ctmp = A[i, 1]*B[j, 1]'
for k = 2:nA
if length(A) == 0 || length(B) == 0
Ctmp = zero(R)
else
Ctmp = zero(A[i, 1]*B[j, 1] + A[i, 1]*B[j, 1])
end
for k = 1:nA
Ctmp += A[i, k]*B[j, k]'
end
C[i,j] = Ctmp
Expand All @@ -519,24 +545,36 @@ function generic_matmatmul!{T,S,R}(C::AbstractVecOrMat{R}, tA, tB, A::AbstractVe
elseif tA == 'T'
if tB == 'N'
for i = 1:mA, j = 1:nB
Ctmp = A[1, i].'B[1, j]
for k = 2:nA
if length(A) == 0 || length(B) == 0
Ctmp = zero(R)
else
Ctmp = zero(A[1, i]*B[1, j] + A[1, i]*B[1, j])
end
for k = 1:nA
Ctmp += A[k, i].'B[k, j]
end
C[i,j] = Ctmp
end
elseif tB == 'T'
for i = 1:mA, j = 1:nB
Ctmp = A[1, i].'B[j, 1].'
for k = 2:nA
if length(A) == 0 || length(B) == 0
Ctmp = zero(R)
else
Ctmp = zero(A[1, i]*B[j, 1] + A[1, i]*B[j, 1])
end
for k = 1:nA
Ctmp += A[k, i].'B[j, k].'
end
C[i,j] = Ctmp
end
else
for i = 1:mA, j = 1:nB
Ctmp = A[1, i].'B[j, 1]'
for k = 2:nA
if length(A) == 0 || length(B) == 0
Ctmp = zero(R)
else
Ctmp = zero(A[1, i]*B[j, 1] + A[1, i]*B[j, 1])
end
for k = 1:nA
Ctmp += A[k, i].'B[j, k]'
end
C[i,j] = Ctmp
Expand All @@ -545,24 +583,36 @@ function generic_matmatmul!{T,S,R}(C::AbstractVecOrMat{R}, tA, tB, A::AbstractVe
else
if tB == 'N'
for i = 1:mA, j = 1:nB
Ctmp = A[1, i]'B[1, j]
for k = 2:nA
if length(A) == 0 || length(B) == 0
Ctmp = zero(R)
else
Ctmp = zero(A[1, i]*B[1, j] + A[1, i]*B[1, j])
end
for k = 1:nA
Ctmp += A[k, i]'B[k, j]
end
C[i,j] = Ctmp
end
elseif tB == 'T'
for i = 1:mA, j = 1:nB
Ctmp = A[1, i]'B[j, 1].'
for k = 2:nA
if length(A) == 0 || length(B) == 0
Ctmp = zero(R)
else
Ctmp = zero(A[1, i]*B[j, 1] + A[1, i]*B[j, 1])
end
for k = 1:nA
Ctmp += A[k, i]'B[j, k].'
end
C[i,j] = Ctmp
end
else
for i = 1:mA, j = 1:nB
Ctmp = A[1, i]'B[j, 1]'
for k = 2:nA
if length(A) == 0 || length(B) == 0
Ctmp = zero(R)
else
Ctmp = zero(A[1, i]*B[j, 1] + A[1, i]*B[j, 1])
end
for k = 1:nA
Ctmp += A[k, i]'B[j, k]'
end
C[i,j] = Ctmp
Expand Down
11 changes: 11 additions & 0 deletions test/linalg/matmul.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,3 +132,14 @@ let A = [1+2im 3+4im; 5+6im 7+8im], B = [2+7im 4+1im; 3+8im 6+5im]
end
end
end

# Issue 11978
A = Array(Matrix{Float64}, 2, 2)
A[1,1] = eye(3)
A[1,2] = eye(3,2)
A[2,1] = eye(2,3)
A[2,2] = eye(2)
b = Array(Vector{Float64}, 2)
b[1] = ones(3)
b[2] = ones(2)
@test A*b == Vector{Float64}[[2,2,1], [2,2]]

0 comments on commit 351e8af

Please sign in to comment.