Skip to content

Commit a77697d

Browse files
DanielVandHDanielVandH
authored
Add * for more matrices defined from layouts (#241)
* Add * for more layouts * bump * cleanup * LTS doesn't have Xoshiro * Missed the infmats * Ambiguities on LTS
1 parent d315ce3 commit a77697d

File tree

5 files changed

+121
-5
lines changed

5 files changed

+121
-5
lines changed

‎Project.toml

+1-1
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "ArrayLayouts"
22
uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
33
authors = ["Sheehan Olver <solver@mac.com>"]
4-
version = "1.10.1"
4+
version = "1.10.2"
55

66
[deps]
77
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

‎src/ArrayLayouts.jl

+2
Original file line numberDiff line numberDiff line change
@@ -219,6 +219,8 @@ end
219219
*(A::Diagonal{<:Any,<:LayoutVector}, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B)
220220
*(A::Diagonal{<:Any,<:LayoutVector}, B::AbstractMatrix) = mul(A, B)
221221
*(A::AbstractMatrix, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B)
222+
*(A::Union{Bidiagonal, Tridiagonal}, B::Diagonal{<:Any, <:LayoutVector}) = mul(A, B) # ambiguity
223+
*(A::Diagonal{<:Any, <:LayoutVector}, B::Union{Bidiagonal, Tridiagonal}) = mul(A, B) # ambiguity
222224
*(A::Diagonal{<:Any,<:LayoutVector}, B::LayoutMatrix) = mul(A, B)
223225
*(A::LayoutMatrix, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B)
224226
*(A::UpperTriangular, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B)

‎src/mul.jl

+19
Original file line numberDiff line numberDiff line change
@@ -362,6 +362,25 @@ end
362362
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::UpperOrLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}) = mul(A, B)
363363
*(A::UpperOrLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}, B::UpperOrLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}) = mul(A, B)
364364

365+
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
366+
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
367+
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
368+
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
369+
*(A::SymTridiagonal, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B) # ambiguity
370+
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
371+
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
372+
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
373+
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
374+
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
375+
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
376+
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
377+
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
378+
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
379+
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
380+
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
381+
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B) # ambiguity
382+
*(A::Diagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
383+
*(A::Diagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
365384

366385
# mul! for subarray of layout matrix
367386
const SubLayoutMatrix = Union{SubArray{<:Any,2,<:LayoutMatrix}, SubArray{<:Any,2,<:AdjOrTrans{<:Any,<:LayoutMatrix}}}

‎test/infinitearrays.jl

+76-2
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,17 @@
11
# Infinite Arrays implementation from
22
# https://github.com/JuliaLang/julia/blob/master/test/testhelpers/InfiniteArrays.jl
33
module InfiniteArrays
4-
using Infinities
5-
export OneToInf
4+
using Infinities, LinearAlgebra, Random
5+
using ..ArrayLayouts: ArrayLayouts, LayoutVector, LayoutMatrix, Mul, DenseColumnMajor
6+
export OneToInf,
7+
InfSymTridiagonal,
8+
InfTridiagonal,
9+
InfBidiagonal,
10+
InfUnitUpperTriangular,
11+
InfUnitLowerTriangular,
12+
InfUpperTriangular,
13+
InfLowerTriangular,
14+
InfDiagonal
615

716
abstract type AbstractInfUnitRange{T<:Real} <: AbstractUnitRange{T} end
817
Base.length(r::AbstractInfUnitRange) = ℵ₀
@@ -48,4 +57,69 @@ module InfiniteArrays
4857
@boundscheck checkbounds(v, first(i))
4958
v[first(i)]:ℵ₀
5059
end
60+
61+
## Methods for testing infinite arrays
62+
struct InfVec{RNG} <: LayoutVector{Float64} # show is broken for InfVec
63+
rng::RNG
64+
data::Vector{Float64}
65+
end
66+
InfVec() = InfVec(copy(Random.seed!(Random.default_rng(), rand(UInt64))), Float64[])
67+
function resizedata!(v::InfVec, i)
68+
n = length(v.data)
69+
i ≤ n && return v
70+
resize!(v.data, i)
71+
for j in (n+1):i
72+
v[j] = rand(v.rng)
73+
end
74+
return v
75+
end
76+
Base.getindex(v::InfVec, i::Int) = (resizedata!(v, i); v.data[i])
77+
Base.setindex!(v::InfVec, r, i::Int) = setindex!(v.data, r, i)
78+
Base.size(v::InfVec) = (ℵ₀,)
79+
Base.axes(v::InfVec) = (OneToInf(),)
80+
ArrayLayouts.MemoryLayout(::Type{<:InfVec}) = DenseColumnMajor()
81+
Base.similar(v::InfVec, ::Type{T}, ::Tuple{OneToInf{Int}}) where {T} = InfVec()
82+
Base.copy(v::InfVec) = InfVec(copy(v.rng), copy(v.data))
83+
84+
struct InfMat{RNG} <: LayoutMatrix{Float64} # show is broken for InfMat
85+
vec::InfVec{RNG}
86+
end
87+
InfMat() = InfMat(InfVec())
88+
function diagtrav_idx(i, j)
89+
band = i + j - 1
90+
nelm = (band * (band - 1)) ÷ 2
91+
return nelm + i
92+
end
93+
Base.getindex(A::InfMat, i::Int, j::Int) = A.vec[diagtrav_idx(i, j)]
94+
Base.setindex!(A::InfMat, r, i::Int, j::Int) = setindex!(A.vec, r, diagtrav_idx(i, j))
95+
Base.size(A::InfMat) = (ℵ₀, ℵ₀)
96+
Base.axes(v::InfMat) = (OneToInf(), OneToInf())
97+
ArrayLayouts.MemoryLayout(::Type{<:InfMat}) = DenseColumnMajor()
98+
Base.copy(A::InfMat) = InfMat(copy(A.vec))
99+
100+
const InfSymTridiagonal = SymTridiagonal{Float64,<:InfVec}
101+
const InfTridiagonal = Tridiagonal{Float64,<:InfVec}
102+
const InfBidiagonal = Bidiagonal{Float64,<:InfVec}
103+
const InfUnitUpperTriangular = UnitUpperTriangular{Float64,<:InfMat}
104+
const InfUnitLowerTriangular = UnitLowerTriangular{Float64,<:InfMat}
105+
const InfUpperTriangular = UpperTriangular{Float64,<:InfMat}
106+
const InfLowerTriangular = LowerTriangular{Float64,<:InfMat}
107+
const InfDiagonal = Diagonal{Float64,<:InfVec}
108+
InfSymTridiagonal() = SymTridiagonal(InfVec(), InfVec())
109+
InfTridiagonal() = Tridiagonal(InfVec(), InfVec(), InfVec())
110+
InfBidiagonal(uplo) = Bidiagonal(InfVec(), InfVec(), uplo)
111+
InfUnitUpperTriangular() = UnitUpperTriangular(InfMat())
112+
InfUnitLowerTriangular() = UnitLowerTriangular(InfMat())
113+
InfUpperTriangular() = UpperTriangular(InfMat())
114+
InfLowerTriangular() = LowerTriangular(InfMat())
115+
InfDiagonal() = Diagonal(InfVec())
116+
Base.copy(D::InfDiagonal) = Diagonal(copy(D.diag))
117+
118+
# 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)
119+
const InfNamedMatrix = Union{InfSymTridiagonal,InfTridiagonal,InfBidiagonal,
120+
InfUnitUpperTriangular,InfUnitLowerTriangular,
121+
InfUpperTriangular,InfLowerTriangular,
122+
InfDiagonal}
123+
const InfMul{L1,L2} = Mul{L1,L2,<:InfNamedMatrix,<:InfNamedMatrix}
124+
Base.copy(M::InfMul{L1,L2}) where {L1,L2} = Mul{L1,L2}(copy(M.A), copy(M.B))
51125
end

‎test/test_layoutarray.jl

+23-2
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
module TestLayoutArray
22

3-
using ArrayLayouts, LinearAlgebra, FillArrays, Test, SparseArrays
4-
using ArrayLayouts: sub_materialize, MemoryLayout, ColumnNorm, RowMaximum, CRowMaximum, @_layoutlmul
3+
using ArrayLayouts, LinearAlgebra, FillArrays, Test, SparseArrays, Random
4+
using ArrayLayouts: sub_materialize, MemoryLayout, ColumnNorm, RowMaximum, CRowMaximum, @_layoutlmul, Mul
55
import ArrayLayouts: triangulardata
66

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

620+
# Tests needed for InfiniteRandomArrays.jl (see https://github.com/DanielVandH/InfiniteRandomArrays.jl/issues/5)
621+
include("infinitearrays.jl")
622+
using .InfiniteArrays
623+
624+
@testset "* for infinite layouts" begin
625+
tup = InfSymTridiagonal(), InfTridiagonal(), InfBidiagonal('U'),
626+
InfBidiagonal('L'),
627+
InfUnitUpperTriangular(), InfUnitLowerTriangular(),
628+
InfUpperTriangular(), InfLowerTriangular(), InfDiagonal();
629+
for (i, A) in enumerate(tup)
630+
A_up, A_lo = A isa Union{UpperTriangular,UnitUpperTriangular}, A isa Union{LowerTriangular,UnitLowerTriangular}
631+
for (j, B) in enumerate(tup)
632+
B_up, B_lo = B isa Union{UpperTriangular,UnitUpperTriangular}, B isa Union{LowerTriangular,UnitLowerTriangular}
633+
((A_up && B_lo) || (A_lo && B_up)) && continue
634+
C = A * B
635+
_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
636+
@test _C ≈ Matrix(A[1:100, 1:102]) * Matrix(B[1:102, 1:100])
637+
end
638+
end
639+
end
640+
620641
end

0 commit comments

Comments
 (0)