Skip to content

Commit 348f5d8

Browse files
jishnubjishnub
authored
Merge branch 'master' into jishnub/diag_mul_bitridiag
2 parents 0e6a009 + 40dc4f3 commit 348f5d8

File tree

5 files changed

+135
-5
lines changed

5 files changed

+135
-5
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
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 <[email protected]>"]
4-
version = "1.10.2"
4+
version = "1.10.4"
55

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

‎src/ArrayLayouts.jl‎

Lines changed: 2 additions & 0 deletions
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‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -347,6 +347,11 @@ end
347347
*(A::Transpose{<:Any,<:LayoutVector}, B::Adjoint{<:Any,<:LayoutMatrix}) = mul(A,B)
348348
*(A::Transpose{<:Any,<:LayoutVector}, B::Transpose{<:Any,<:LayoutMatrix}) = mul(A,B)
349349

350+
# Disambiguation with FillArrays
351+
*(A::AbstractFill{<:Any,2}, B::LayoutVector) = invoke(*, Tuple{AbstractFill{<:Any,2}, AbstractVector}, A, B)
352+
*(A::Adjoint{<:Any, <:LayoutVector}, B::AbstractFill{<:Any,2}) = invoke(*, Tuple{Adjoint{<:Any, <:AbstractVector}, AbstractFill{<:Any,2}}, A, B)
353+
*(A::Transpose{<:Any, <:LayoutVector}, B::AbstractFill{<:Any,2}) = invoke(*, Tuple{Transpose{<:Any, <:AbstractVector}, AbstractFill{<:Any,2}}, A, B)
354+
350355
## special routines introduced in v0.9. We need to avoid these to support ∞-arrays
351356

352357
*(x::Adjoint{<:Any,<:LayoutVector}, D::Diagonal{<:Any,<:LayoutVector}) = mul(x, D)
@@ -362,6 +367,25 @@ end
362367
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::UpperOrLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}) = mul(A, B)
363368
*(A::UpperOrLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}, B::UpperOrLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}) = mul(A, B)
364369

370+
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
371+
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
372+
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
373+
*(A::SymTridiagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
374+
*(A::SymTridiagonal, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B) # ambiguity
375+
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
376+
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
377+
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
378+
*(A::Tridiagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
379+
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
380+
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
381+
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
382+
*(A::Bidiagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
383+
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
384+
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::Tridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
385+
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::Bidiagonal{<:Any,<:LayoutVector}) = mul(A, B)
386+
*(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B) # ambiguity
387+
*(A::Diagonal{<:Any,<:LayoutVector}, B::SymTridiagonal{<:Any,<:LayoutVector}) = mul(A, B)
388+
*(A::Diagonal{<:Any,<:LayoutVector}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B)
365389

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

‎test/infinitearrays.jl‎

Lines changed: 76 additions & 2 deletions
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‎

Lines changed: 32 additions & 2 deletions
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
import LinearAlgebra: Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal
77

@@ -642,4 +642,34 @@ triangulardata(A::MyUpperTriangular) = triangulardata(A.A)
642642
VERSION >= v"1.9" && @test U / MyMatrix(A) ≈ U / A
643643
end
644644

645+
# Tests needed for InfiniteRandomArrays.jl (see https://github.com/DanielVandH/InfiniteRandomArrays.jl/issues/5)
646+
include("infinitearrays.jl")
647+
using .InfiniteArrays
648+
649+
@testset "* for infinite layouts" begin
650+
tup = InfSymTridiagonal(), InfTridiagonal(), InfBidiagonal('U'),
651+
InfBidiagonal('L'),
652+
InfUnitUpperTriangular(), InfUnitLowerTriangular(),
653+
InfUpperTriangular(), InfLowerTriangular(), InfDiagonal();
654+
for (i, A) in enumerate(tup)
655+
A_up, A_lo = A isa Union{UpperTriangular,UnitUpperTriangular}, A isa Union{LowerTriangular,UnitLowerTriangular}
656+
for (j, B) in enumerate(tup)
657+
B_up, B_lo = B isa Union{UpperTriangular,UnitUpperTriangular}, B isa Union{LowerTriangular,UnitLowerTriangular}
658+
((A_up && B_lo) || (A_lo && B_up)) && continue
659+
C = A * B
660+
_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
661+
@test _C ≈ Matrix(A[1:100, 1:102]) * Matrix(B[1:102, 1:100])
662+
end
663+
end
664+
end
665+
666+
@testset "disambiguation with FillArrays" begin
667+
v = [1,2,3]
668+
lv = MyVector(v)
669+
F = Fill(2, 3, 3)
670+
@test F * lv == F * v
671+
@test lv' * F == v' * F
672+
@test transpose(lv) * F == transpose(v) * F
673+
end
674+
645675
end

0 commit comments

Comments
 (0)