*(::AbstractArray, ::AbstractArray) is wrong

[odow]

x-ref #174

These definitions place some very strong assumptions on the input types that are not present in *(A, B):

MutableArithmetics.jl/src/implementations/LinearAlgebra.jl

Lines 380 to 403 in 201aafb

	# Does what `LinearAlgebra/src/matmul.jl` does for abstract matrices and
	# vectors: estimate the resulting element type, allocate the resulting array but
	# it redirects to `mul_to!` instead of `LinearAlgebra.mul!`.
	function operate(
	::typeof(*),
	A::AbstractMatrix{S},
	B::AbstractVector{T},
	) where {T,S}
	C = undef_array(promote_array_mul(typeof(A), typeof(B)), axes(A, 1))
	return operate_to!(C, *, A, B)
	end
	function operate(
	::typeof(*),
	A::AbstractMatrix{S},
	B::AbstractMatrix{T},
	) where {T,S}
	C = undef_array(
	promote_array_mul(typeof(A), typeof(B)),
	axes(A, 1),
	axes(B, 2),
	)
	return operate_to!(C, *, A, B)
	end

For example, none of these work:

@testset "Abstract eltype in matmul" begin
    for T in (Any, Union{String,Int})
        x, x12, x22 = T[1, 2], T[1 2], T[1 2; 3 4]
        @test MA.operate(*, x, x') == x * x'
        @test MA.operate(*, x', x) == x' * x
        @test MA.operate(*, x12, x) == x12 * x
        @test MA.operate(*, x22, x) == x22 * x
        @test MA.operate(*, x', x22) == x' * x22
        @test MA.operate(*, x12, x22) == x12 * x22
        @test MA.operate(*, x22, x22) == x22 * x22
    end
end

[blegat]

We could replace typeof with _concrete_typeof similar to _concrete_eltype of #171

[odow]

I still think that there's an issue we can't resolve:

julia> y = Vector{Union{Int,Float64}}([1, 1.5])
2-element Vector{Union{Float64, Int64}}:
 1.0
 1.5

julia> y = Union{Int,Float64}[1, 1.5]
2-element Vector{Union{Float64, Int64}}:
 1
 1.5

julia> y * y'
2×2 Matrix{Real}:
 1    1.5
 1.5  2.25

julia> y' * y
3.25

Base doesn't rely on the container types at compile time. It uses promotion based on the run-time values. Pre-allocating C cannot work in all cases that we expect it to, unless we have an additional loop over all of A and B?

[odow]

I have a very simple fix that I'm just testing 😄

[odow]

Urgh. Even Base can't get this right:

julia> M = LinearAlgebra.LowerTriangular
LowerTriangular

julia> T = Union{Int,Float64}
Union{Float64, Int64}

julia> x = M(T[1 2.5; 3.5 4])
2×2 LowerTriangular{Union{Float64, Int64}, Matrix{Union{Float64, Int64}}}:
 1     ⋅ 
 3.5  4

julia> x * x
ERROR: InexactError: Int64(3.5)
Stacktrace:
 [1] Int64
   @ ./float.jl:723 [inlined]
 [2] convert
   @ ./number.jl:7 [inlined]
 [3] setindex!
   @ ./array.jl:843 [inlined]
 [4] copyto_unaliased!(deststyle::IndexLinear, dest::Matrix{Int64}, srcstyle::IndexCartesian, src::LowerTriangular{Union{Float64, Int64}, Matrix{Union{Float64, Int64}}})
   @ Base ./abstractarray.jl:976
 [5] copyto!
   @ ./abstractarray.jl:950 [inlined]
 [6] *(A::LowerTriangular{Union{Float64, Int64}, Matrix{Union{Float64, Int64}}}, B::LowerTriangular{Union{Float64, Int64}, Matrix{Union{Float64, Int64}}})
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:1522
 [7] top-level scope
   @ REPL[17]:1

[blegat]

Base doesn't rely on the container types at compile time. It uses promotion based on the run-time values. Pre-allocating C cannot work in all cases that we expect it to, unless we have an additional loop over all of A and B?

Yes, LinearAlgebra will start creating a Matrix{Int} and then when floating points arrive, it will broaden the eltype so that the Int still fit instead of promoting and converting the Ints already produced into Float64 (since that might not always work).