Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add * for more matrices defined from layouts #241

Merged
merged 6 commits into from
Jul 7, 2024
Merged

Add * for more matrices defined from layouts #241

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
4304f5e
Add * for more layouts
DanielVandH Jul 5, 2024
b7c2b8c
bump
DanielVandH Jul 5, 2024
88abe00
cleanup
DanielVandH Jul 5, 2024
27f1446
LTS doesn't have Xoshiro
DanielVandH Jul 5, 2024
b04f5bf
Missed the infmats
DanielVandH Jul 5, 2024
33a160e
Ambiguities on LTS
DanielVandH Jul 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "ArrayLayouts"
uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
authors = ["Sheehan Olver <solver@mac.com>"]
version = "1.10.1"
version = "1.10.2"

[deps]
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
Expand Down
2 changes: 2 additions & 0 deletions src/ArrayLayouts.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -219,6 +219,8 @@ end
*(A::Diagonal{<:Any,<:LayoutVector}, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::Diagonal{<:Any,<:LayoutVector}, B::AbstractMatrix) = mul(A, B)
*(A::AbstractMatrix, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::Union{Bidiagonal, Tridiagonal}, B::Diagonal{<:Any, <:LayoutVector}) = mul(A, B) # ambiguity
*(A::Diagonal{<:Any, <:LayoutVector}, B::Union{Bidiagonal, Tridiagonal}) = mul(A, B) # ambiguity
*(A::Diagonal{<:Any,<:LayoutVector}, B::LayoutMatrix) = mul(A, B)
*(A::LayoutMatrix, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::UpperTriangular, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B)
Expand Down
19 changes: 19 additions & 0 deletions src/mul.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -362,6 +362,25 @@ end
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::UpperOrLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}) = mul(A, B)
*(A::UpperOrLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}, B::UpperOrLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}) = mul(A, B)

*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
*(A::SymTridiagonal, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B) # ambiguity
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B) # ambiguity
*(A::Diagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
*(A::Diagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
Copy link
Member

@jishnub jishnub Jul 5, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

@jishnub jishnub Jul 5, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


# mul! for subarray of layout matrix
const SubLayoutMatrix = Union{SubArray{<:Any,2,<:LayoutMatrix}, SubArray{<:Any,2,<:AdjOrTrans{<:Any,<:LayoutMatrix}}}
Expand Down
78 changes: 76 additions & 2 deletions test/infinitearrays.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,17 @@
# Infinite Arrays implementation from
# https://github.com/JuliaLang/julia/blob/master/test/testhelpers/InfiniteArrays.jl
module InfiniteArrays
using Infinities
export OneToInf
using Infinities, LinearAlgebra, Random
using ..ArrayLayouts: ArrayLayouts, LayoutVector, LayoutMatrix, Mul, DenseColumnMajor
export OneToInf,
InfSymTridiagonal,
InfTridiagonal,
InfBidiagonal,
InfUnitUpperTriangular,
InfUnitLowerTriangular,
InfUpperTriangular,
InfLowerTriangular,
InfDiagonal

abstract type AbstractInfUnitRange{T<:Real} <: AbstractUnitRange{T} end
Base.length(r::AbstractInfUnitRange) = ℵ₀
Expand Down Expand Up @@ -48,4 +57,69 @@ module InfiniteArrays
@boundscheck checkbounds(v, first(i))
v[first(i)]:ℵ₀
end

## Methods for testing infinite arrays
struct InfVec{RNG} <: LayoutVector{Float64} # show is broken for InfVec
rng::RNG
data::Vector{Float64}
end
InfVec() = InfVec(copy(Random.seed!(Random.default_rng(), rand(UInt64))), Float64[])
function resizedata!(v::InfVec, i)
n = length(v.data)
i ≤ n && return v
resize!(v.data, i)
for j in (n+1):i
v[j] = rand(v.rng)
end
return v
end
Base.getindex(v::InfVec, i::Int) = (resizedata!(v, i); v.data[i])
Base.setindex!(v::InfVec, r, i::Int) = setindex!(v.data, r, i)
Base.size(v::InfVec) = (ℵ₀,)
Base.axes(v::InfVec) = (OneToInf(),)
ArrayLayouts.MemoryLayout(::Type{<:InfVec}) = DenseColumnMajor()
Base.similar(v::InfVec, ::Type{T}, ::Tuple{OneToInf{Int}}) where {T} = InfVec()
Base.copy(v::InfVec) = InfVec(copy(v.rng), copy(v.data))

struct InfMat{RNG} <: LayoutMatrix{Float64} # show is broken for InfMat
vec::InfVec{RNG}
end
InfMat() = InfMat(InfVec())
function diagtrav_idx(i, j)
band = i + j - 1
nelm = (band * (band - 1)) ÷ 2
return nelm + i
end
Base.getindex(A::InfMat, i::Int, j::Int) = A.vec[diagtrav_idx(i, j)]
Base.setindex!(A::InfMat, r, i::Int, j::Int) = setindex!(A.vec, r, diagtrav_idx(i, j))
Base.size(A::InfMat) = (ℵ₀, ℵ₀)
Base.axes(v::InfMat) = (OneToInf(), OneToInf())
ArrayLayouts.MemoryLayout(::Type{<:InfMat}) = DenseColumnMajor()
Base.copy(A::InfMat) = InfMat(copy(A.vec))

const InfSymTridiagonal = SymTridiagonal{Float64,<:InfVec}
const InfTridiagonal = Tridiagonal{Float64,<:InfVec}
const InfBidiagonal = Bidiagonal{Float64,<:InfVec}
const InfUnitUpperTriangular = UnitUpperTriangular{Float64,<:InfMat}
const InfUnitLowerTriangular = UnitLowerTriangular{Float64,<:InfMat}
const InfUpperTriangular = UpperTriangular{Float64,<:InfMat}
const InfLowerTriangular = LowerTriangular{Float64,<:InfMat}
const InfDiagonal = Diagonal{Float64,<:InfVec}
InfSymTridiagonal() = SymTridiagonal(InfVec(), InfVec())
InfTridiagonal() = Tridiagonal(InfVec(), InfVec(), InfVec())
InfBidiagonal(uplo) = Bidiagonal(InfVec(), InfVec(), uplo)
InfUnitUpperTriangular() = UnitUpperTriangular(InfMat())
InfUnitLowerTriangular() = UnitLowerTriangular(InfMat())
InfUpperTriangular() = UpperTriangular(InfMat())
InfLowerTriangular() = LowerTriangular(InfMat())
InfDiagonal() = Diagonal(InfVec())
Base.copy(D::InfDiagonal) = Diagonal(copy(D.diag))

# Without LazyArrays we have no access to the lazy machinery, so we must define copy(::Mul) to leave mul(A, B) as a lazy Mul(A, B)
const InfNamedMatrix = Union{InfSymTridiagonal,InfTridiagonal,InfBidiagonal,
InfUnitUpperTriangular,InfUnitLowerTriangular,
InfUpperTriangular,InfLowerTriangular,
InfDiagonal}
const InfMul{L1,L2} = Mul{L1,L2,<:InfNamedMatrix,<:InfNamedMatrix}
Base.copy(M::InfMul{L1,L2}) where {L1,L2} = Mul{L1,L2}(copy(M.A), copy(M.B))
end
25 changes: 23 additions & 2 deletions test/test_layoutarray.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
module TestLayoutArray

using ArrayLayouts, LinearAlgebra, FillArrays, Test, SparseArrays
using ArrayLayouts: sub_materialize, MemoryLayout, ColumnNorm, RowMaximum, CRowMaximum, @_layoutlmul
using ArrayLayouts, LinearAlgebra, FillArrays, Test, SparseArrays, Random
using ArrayLayouts: sub_materialize, MemoryLayout, ColumnNorm, RowMaximum, CRowMaximum, @_layoutlmul, Mul
import ArrayLayouts: triangulardata

struct MyMatrix <: LayoutMatrix{Float64}
Expand Down Expand Up @@ -617,4 +617,25 @@ triangulardata(A::MyUpperTriangular) = triangulardata(A.A)
VERSION >= v"1.9" && @test U / MyMatrix(A) ≈ U / A
end

# Tests needed for InfiniteRandomArrays.jl (see https://github.com/DanielVandH/InfiniteRandomArrays.jl/issues/5)
include("infinitearrays.jl")
using .InfiniteArrays

@testset "* for infinite layouts" begin
tup = InfSymTridiagonal(), InfTridiagonal(), InfBidiagonal('U'),
InfBidiagonal('L'),
InfUnitUpperTriangular(), InfUnitLowerTriangular(),
InfUpperTriangular(), InfLowerTriangular(), InfDiagonal();
for (i, A) in enumerate(tup)
A_up, A_lo = A isa Union{UpperTriangular,UnitUpperTriangular}, A isa Union{LowerTriangular,UnitLowerTriangular}
for (j, B) in enumerate(tup)
B_up, B_lo = B isa Union{UpperTriangular,UnitUpperTriangular}, B isa Union{LowerTriangular,UnitLowerTriangular}
((A_up && B_lo) || (A_lo && B_up)) && continue
C = A * B
_C = [C[i, j] for i in 1:100, j in 1:100] # else we need to fix the C[1:100, 1:100] MethodError from _getindex(::Mul, ...). This is easier
@test _C ≈ Matrix(A[1:100, 1:102]) * Matrix(B[1:102, 1:100])
end
end
end

end
Loading