Add * for more matrices defined from layouts

[DanielVandH]

Need this to support multiplication between infinite matrices, as found in DanielVandH/InfiniteRandomArrays.jl#5.

To make sure the methods are being tested properly, I extended the minimal InfiniteArrays.jl module in the tests to include a (hopefully minimal) definition of an infinite vector and an infinite matrix, and the resulting named matrices to mimic those from InfiniteRandomMatrices.jl.

InfVec and InfMat were initially just matrices which returned the indices, e.g. InfVec()[i] = i, but I wanted to be sure the methods were properly copying the data and accessing the correct entries, so I convert it into a random vector.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 90.53%. Comparing base (d315ce3) to head (33a160e).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #241      +/-   ##
==========================================
+ Coverage   90.43%   90.53%   +0.10%     
==========================================
  Files          11       11              
  Lines        1902     1923      +21     
==========================================
+ Hits         1720     1741      +21     
  Misses        182      182              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[DanielVandH]

I'll deal with LTS later

[jishnub]

This should forward to *(A::Diagonal, B.data::LayoutMatrix), so do we need to specialize the method for the triangular wrapper? Or will specializing it for LayoutMatrix suffice? *(::Diagonal, ::LayoutMatrix) method already exists, which should also make this work?

[DanielVandH]

Are you referring to all four of those methods? I think actually the second and fourth methods in your highlight aren't needed, so I could remove those. Otherwise I'm not sure what you mean

[DanielVandH]

Actually, they are all needed due to ambiguities that these fixes have to introduce

[jishnub]

I'm referring to the *(A::Diagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) method on line 383, and the other order. Does this lead to ambiguities? If it does, I think this should be resolved in LinearAlgebra as well.

[DanielVandH]

This is the complete error from the ambiguity that made me add that line (the other direction doesn't have an ambiguity)

julia> InfDiagonal() * InfUnitUpperTriangular() # has to be unit to get the ambiguity
ERROR: MethodError: *(::Diagonal{Float64, Main.InfiniteArrays.InfVec{Xoshiro}}, ::UnitUpperTriangular{Float64, Main.InfiniteArrays.InfMat{Xoshiro}}) is ambiguous.

Candidates:
  *(D::Diagonal, A::UnitUpperTriangular)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\diagonal.jl:578       
  *(A::Diagonal{<:Any, <:LayoutVector}, B::AbstractMatrix)
    @ ArrayLayouts c:\Users\djv23\.julia\dev\ArrayLayouts.jl\src\ArrayLayouts.jl:220
  *(D::Diagonal, A::LinearAlgebra.AbstractTriangular)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\diagonal.jl:835       
  *(A::AbstractMatrix, B::LinearAlgebra.AbstractTriangular)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:1500    
  *(D::Diagonal, A::AbstractMatrix)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\diagonal.jl:292       

Possible fix, define
  *(::Diagonal{T, <:LayoutVector{T}} where T, ::UnitUpperTriangular)

Stacktrace:
 [1] top-level scope
   @ REPL[16]:1

I could make the method be Union{UnitUpperTriangular, UnitLowerTriangular} instead.

[DanielVandH]

I could make the method be Union{UnitUpperTriangular, UnitLowerTriangular} instead.

Actually nevermind, that just adds another ambiguity that forces me to add the more general line.. now I remember.

[jishnub]

Could you provide a MWE? How are InfUnitUpperTriangular and InfDiagonal defined? Ideally, LinearAlgebra should not be introducing these ambiguities, but we may have to add these methods for now, and version-limit these later.

[DanielVandH]

They are defined in this PR, built from my InfVec and InfMat definitions. Without this PR, you can reproduce it using this MWE:

using ArrayLayouts, LinearAlgebra
struct MyVec <: LayoutVector{Float64} 
    A::Vector{Float64}
end
struct MyMat <: LayoutMatrix{Float64}
    A::Matrix{Float64}
end
Base.size(v::MyVec) = size(v.A)
Base.getindex(v::MyVec, i::Int) = v.A[i] 
Base.size(v::MyMat) = size(v.A)
Base.getindex(v::MyMat, i::Int, j::Int) = v.A[i, j]
D = Diagonal(MyVec(rand(10)))
A = UnitUpperTriangular(MyMat(rand(10, 10)))
D * A

julia> D * A
ERROR: MethodError: *(::Diagonal{Float64, MyVec}, ::UnitUpperTriangular{Float64, MyMat}) is ambiguous.

Candidates:
  *(D::Diagonal, A::UnitUpperTriangular)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\diagonal.jl:578
  *(A::Diagonal{<:Any, <:LayoutVector}, B::AbstractMatrix)
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\VzDAX\src\ArrayLayouts.jl:220
  *(D::Diagonal, A::LinearAlgebra.AbstractTriangular)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\diagonal.jl:835
  *(A::AbstractMatrix, B::LinearAlgebra.AbstractTriangular)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:1500
  *(D::Diagonal, A::AbstractMatrix)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\diagonal.jl:292

Possible fix, define
  *(::Diagonal{T, <:LayoutVector{T}} where T, ::UnitUpperTriangular)

Stacktrace:
 [1] top-level scope
   @ Untitled-1:14

[DanielVandH]

Is there anything else you want for this to be merged @jishnub? I'm not sure what would be defined on LinearAlgebra.jl's side but I think this would be good to go for now.