From 8fa37c8ffaf51628fe0e2d7815756c30f0cfd758 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 19 Jan 2019 21:28:01 +0000 Subject: [PATCH 01/15] start applied --- src/LazyArrays.jl | 4 +- src/lazyapplying.jl | 20 ++++++++++ src/linalg/blasbroadcasting.jl | 12 +++--- src/linalg/blasmul.jl | 8 ++-- src/linalg/inv.jl | 20 +++++----- src/linalg/mul.jl | 70 ++++++++++++---------------------- test/applytests.jl | 12 ++++++ test/multests.jl | 3 -- test/runtests.jl | 1 + 9 files changed, 80 insertions(+), 70 deletions(-) create mode 100644 src/lazyapplying.jl create mode 100644 test/applytests.jl diff --git a/src/LazyArrays.jl b/src/LazyArrays.jl index 721d62d7..911aaf1b 100644 --- a/src/LazyArrays.jl +++ b/src/LazyArrays.jl @@ -39,10 +39,12 @@ import FillArrays: AbstractFill import StaticArrays: StaticArrayStyle export Mul, MulArray, MulVector, MulMatrix, InvMatrix, PInvMatrix, - Hcat, Vcat, Kron, BroadcastArray, cache, Ldiv, Inv, PInv, Diff, Cumsum + Hcat, Vcat, Kron, BroadcastArray, cache, Ldiv, Inv, PInv, Diff, Cumsum, + applied include("memorylayout.jl") include("cache.jl") +include("lazyapplying.jl") include("lazybroadcasting.jl") include("lazyconcat.jl") include("linalg/linalg.jl") diff --git a/src/lazyapplying.jl b/src/lazyapplying.jl new file mode 100644 index 00000000..4a34a331 --- /dev/null +++ b/src/lazyapplying.jl @@ -0,0 +1,20 @@ + + + + +abstract type ApplyStyle end +struct DefaultApplyStyle <: ApplyStyle end +struct LayoutApplyStyle{Layouts<:Tuple} <: ApplyStyle + layouts::Layouts +end + +ApplyStyle(f, args...) = DefaultApplyStyle() + +struct Applied{Style<:ApplyStyle, F, Args<:Tuple} + style::Style + f::F + args::Args +end + +applied(f, args...) = Applied(ApplyStyle(f, args...), f, args) +materialize(A::Applied{DefaultApplyStyle}) = A.f(materialize.(A.args)...) diff --git a/src/linalg/blasbroadcasting.jl b/src/linalg/blasbroadcasting.jl index 4c33502f..28b46fbc 100644 --- a/src/linalg/blasbroadcasting.jl +++ b/src/linalg/blasbroadcasting.jl @@ -61,7 +61,7 @@ instantiate(bc::Broadcasted{<:ArrayMulArrayStyle}) = bc @inline function _copyto!(_, dest::AbstractArray{T}, M::ArrayMulArray) where T - A,B = M.factors + A,B = M.args materialize!(MulAdd(one(T), A, B, zero(T), dest)) end @@ -73,34 +73,34 @@ end @inline function _copyto!(_, dest::AbstractArray{T}, bc::BConstArrayMulArray) where T α,M = bc.args - A,B = M.factors + A,B = M.args materialize!(MulAdd(α, A, B, zero(T), dest)) end @inline function _copyto!(_, dest::AbstractArray{T}, bc::BArrayMulArrayPlusArray) where T M,C = bc.args - A,B = M.factors + A,B = M.args copyto!(dest, MulAdd(one(T), A, B, one(T), C)) end @inline function _copyto!(_, dest::AbstractArray{T}, bc::BArrayMulArrayPlusConstArray) where T M,βc = bc.args β,C = βc.args - A,B = M.factors + A,B = M.args copyto!(dest, MulAdd(one(T), A, B, β, C)) end @inline function _copyto!(_, dest::AbstractArray{T}, bc::BConstArrayMulArrayPlusArray) where T αM,C = bc.args α,M = αM.args - A,B = M.factors + A,B = M.args copyto!(dest, MulAdd(α, A, B, one(T), C)) end @inline function _copyto!(_, dest::AbstractArray, bc::BConstArrayMulArrayPlusConstArray) αM,βc = bc.args α,M = αM.args - A,B = M.factors + A,B = M.args β,C = βc.args copyto!(dest, MulAdd(α, A, B, β, C)) end diff --git a/src/linalg/blasmul.jl b/src/linalg/blasmul.jl index adadbe0f..79908d9f 100644 --- a/src/linalg/blasmul.jl +++ b/src/linalg/blasmul.jl @@ -329,7 +329,7 @@ materialize!(M::BlasMatMulVec{<:HermitianLayout{<:AbstractRowMajor},<:AbstractSt @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatMulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractColumnMajor}, <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat} - A,x = M.factors + A,x = M.args x ≡ dest || copyto!(dest, x) BLAS.trmv!(UPLO, 'N', UNIT, triangulardata(A), dest) end @@ -337,7 +337,7 @@ end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatMulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractRowMajor}, <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat} - A,x = M.factors + A,x = M.args x ≡ dest || copyto!(dest, x) BLAS.trmv!(UPLO, 'T', UNIT, transpose(triangulardata(A)), dest) end @@ -346,7 +346,7 @@ end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatMulVec{<:TriangularLayout{UPLO,UNIT,<:ConjLayout{<:AbstractRowMajor}}, <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat} - A,x = M.factors + A,x = M.args x ≡ dest || copyto!(dest, x) BLAS.trmv!(UPLO, 'C', UNIT, triangulardata(A)', dest) end @@ -356,7 +356,7 @@ end function _copyto!(_, dest::AbstractMatrix, M::MatMulMat{<:TriangularLayout}) - A,X = M.factors + A,X = M.args size(dest,2) == size(X,2) || thow(DimensionMismatch("Dimensions must match")) @views for j in axes(dest,2) dest[:,j] .= Mul(A, X[:,j]) diff --git a/src/linalg/inv.jl b/src/linalg/inv.jl index fb75169d..7d2d110f 100644 --- a/src/linalg/inv.jl +++ b/src/linalg/inv.jl @@ -85,7 +85,7 @@ macro lazyldiv(Typ) end *(A::AbstractPInv, B, C...) = materialize(Mul(A,B, C...)) -*(A::AbstractPInv, B::Mul) = materialize(Mul(A, B.factors...)) +*(A::AbstractPInv, B::Mul) = materialize(Mul(A, B.args...)) similar(A::AbstractPInv, ::Type{T}) where T = Array{T}(undef, size(A)) similar(M::ArrayLdivArray, ::Type{T}) where T = Array{T}(undef, size(M)) @@ -99,12 +99,12 @@ end if VERSION ≥ v"1.1-pre" function _copyto!(_, dest::AbstractArray, M::ArrayLdivArray) - Ai, B = M.factors + Ai, B = M.args ldiv!(dest, factorize(pinv(Ai)), B) end else function _copyto!(_, dest::AbstractArray, M::ArrayLdivArray) - Ai, B = M.factors + Ai, B = M.args ldiv!(dest, factorize(pinv(Ai)), copy(B)) end end @@ -120,7 +120,7 @@ broadcastable(M::MatLdivVec) = M ### function _copyto!(_, dest::AbstractArray, M::ArrayLdivArray{<:TriangularLayout}) - Ai, B = M.factors + Ai, B = M.args dest ≡ B || (dest .= B) ldiv!(pinv(Ai), dest) end @@ -128,7 +128,7 @@ end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{<:TriangularLayout{UPLO,UNIT,<:AbstractColumnMajor}, <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat} - Ai,B = M.factors + Ai,B = M.args B ≡ dest || copyto!(dest, B) BLAS.trsv!(UPLO, 'N', UNIT, triangulardata(pinv(Ai)), dest) end @@ -136,7 +136,7 @@ end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{<:TriangularLayout{'U',UNIT,<:AbstractRowMajor}, <:AbstractStridedLayout, T, T}) where {UNIT,T <: BlasFloat} - Ai,B = M.factors + Ai,B = M.args B ≡ dest || copyto!(dest, B) BLAS.trsv!('L', 'T', UNIT, transpose(triangulardata(pinv(Ai))), dest) end @@ -144,7 +144,7 @@ end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:AbstractRowMajor}, <:AbstractStridedLayout, T, T}) where {UNIT,T <: BlasFloat} - Ai,B = M.factors + Ai,B = M.args B ≡ dest || copyto!(dest, B) BLAS.trsv!('U', 'T', UNIT, transpose(triangulardata(pinv(Ai))), dest) end @@ -153,7 +153,7 @@ end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{T, <:TriangularLayout{'U',UNIT,<:ConjLayout{<:AbstractRowMajor}}, <:AbstractStridedLayout, T, T}) where {UNIT,T <: BlasFloat} - Ai,B = M.factors + Ai,B = M.args B ≡ dest || copyto!(dest, B) BLAS.trsv!('L', 'C', UNIT, triangulardata(pinv(Ai))', dest) end @@ -161,13 +161,13 @@ end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:ConjLayout{<:AbstractRowMajor}}, <:AbstractStridedLayout, T, T}) where {UNIT,T <: BlasFloat} - Ai,B = M.factors + Ai,B = M.args B ≡ dest || copyto!(dest, B) BLAS.trsv!('U', 'C', UNIT, triangulardata(pinv(Ai))', dest) end function _copyto!(_, dest::AbstractMatrix, M::MatLdivMat{<:TriangularLayout}) - A,X = M.factors + A,X = M.args size(dest,2) == size(X,2) || thow(DimensionMismatch("Dimensions must match")) @views for j in axes(dest,2) dest[:,j] .= Mul(A, X[:,j]) diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index a69ab341..d89433ef 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -7,42 +7,20 @@ function checkdimensions(A, B, C...) checkdimensions(B, C...) end -struct Mul{Styles<:Tuple, Factors<:Tuple} - styles::Styles - factors::Factors - function Mul{S,F}(styles::S, factors::F) where {S,F} - checkdimensions(factors...) - new{S,F}(styles,factors) - end -end - -Mul(styles::S, factors::F) where {S<:Tuple,F<:Tuple} = Mul{S,F}(styles, factors) - -Mul(A::Tuple) = Mul(MemoryLayout.(A), A) - -""" - Mul(A1, A2, …, AN) - -represents lazy multiplication A1*A2*…*AN. The factors must have compatible axes. -If any argument is itself a Mul, it automatically gets flatten. That is, -we assume associativity. Use Mul((A, B, C)) to stop flattening -""" -Mul(A...) = flatten(Mul(A)) +const Mul{Styles<:Tuple, Factors<:Tuple} = Applied{<:LayoutApplyStyle{Styles}, typeof(*), Factors} -_flatten() = () -_flatten(A, B...) = (A, _flatten(B...)...) -_flatten(A::Mul, B...) = _flatten(A.factors..., B...) -flatten(A::Mul) = Mul(_flatten(A.factors...)) +ApplyStyle(::typeof(*), args::AbstractArray...) = LayoutApplyStyle(MemoryLayout.(args)) +Mul(A...) = applied(*, A...) const Mul2{StyleA, StyleB, AType, BType} = Mul{<:Tuple{StyleA,StyleB}, <:Tuple{AType,BType}} _mul_eltype(a) = eltype(a) _mul_eltype(a, b...) = Base.promote_op(*, eltype(a), _mul_eltype(b...)) -eltype(M::Mul) = _mul_eltype(M.factors...) +eltype(M::Mul) = _mul_eltype(M.args...) size(M::Mul, p::Int) = size(M)[p] axes(M::Mul, p::Int) = axes(M)[p] -ndims(M::Mul) = ndims(last(M.factors)) +ndims(M::Mul) = ndims(last(M.args)) length(M::Mul) = prod(size(M)) size(M::Mul) = length.(axes(M)) @@ -50,7 +28,7 @@ size(M::Mul) = length.(axes(M)) _mul_axes(ax1, ::Tuple{}) = (ax1,) _mul_axes(ax1, ::Tuple{<:Any}) = (ax1,) _mul_axes(ax1, (_,ax2)::Tuple{<:Any,<:Any}) = (ax1,ax2) -axes(M::Mul) = _mul_axes(axes(first(M.factors),1), axes(last(M.factors))) +axes(M::Mul) = _mul_axes(axes(first(M.args),1), axes(last(M.args))) axes(M::Mul{Tuple{}}) = () similar(M::Mul) = similar(M, eltype(M)) @@ -62,25 +40,25 @@ similar(M::Mul) = similar(M, eltype(M)) materializes arrays iteratively, left-to-right. """ -lmaterialize(M::Mul) = _lmaterialize(M.factors...) +lmaterialize(M::Mul) = _lmaterialize(M.args...) -_lmaterialize(A, B) = materialize(Mul((A,B))) -_lmaterialize(A, B, C, D...) = _lmaterialize(materialize(Mul((A,B))), C, D...) +_lmaterialize(A, B) = materialize(Mul(A,B)) +_lmaterialize(A, B, C, D...) = _lmaterialize(materialize(Mul(A,B)), C, D...) """ rmaterialize(M::Mul) materializes arrays iteratively, right-to-left. """ -rmaterialize(M::Mul) = _rmaterialize(reverse(M.factors)...) +rmaterialize(M::Mul) = _rmaterialize(reverse(M.args)...) -_rmaterialize(Z, Y) = materialize(Mul((Y,Z))) -_rmaterialize(Z, Y, X, W...) = _rmaterialize(materialize(Mul((Y,Z))), X, W...) +_rmaterialize(Z, Y) = materialize(Mul(Y,Z)) +_rmaterialize(Z, Y, X, W...) = _rmaterialize(materialize(Mul(Y,Z)), X, W...) -*(A::Mul, B::Mul) = materialize(Mul(A.factors..., B.factors...)) -*(A::Mul, B) = materialize(Mul(A.factors..., B)) -*(A, B::Mul) = materialize(Mul(A, B.factors...)) +# *(A::Mul, B::Mul) = materialize(Mul(A.args..., B.args...)) +# *(A::Mul, B) = materialize(Mul(A.args..., B)) +# *(A, B::Mul) = materialize(Mul(A, B.args...)) ⋆(A...) = Mul(A...) @@ -130,7 +108,7 @@ colsupport(A, j) = colsupport(MemoryLayout(A), A, j) function getindex(M::MatMulVec, k::Integer) - A,B = M.factors + A,B = M.args ret = zero(eltype(M)) for j = rowsupport(A, k) ret += A[k,j] * B[j] @@ -151,7 +129,7 @@ getindex(M::MatMulVec, k::CartesianIndex{1}) = M[convert(Int, k)] const MatMulMat{styleA, styleB, T, V} = ArrayMulArray{styleA, styleB, 2, 2, T, V} function getindex(M::MatMulMat, k::Integer, j::Integer) - A,B = M.factors + A,B = M.args ret = zero(eltype(M)) @inbounds for ℓ in (rowsupport(A,k) ∩ colsupport(B,j)) ret += A[k,ℓ] * B[ℓ,j] @@ -167,7 +145,7 @@ getindex(M::MatMulMat, kj::CartesianIndex{2}) = M[kj[1], kj[2]] ##### function getindex(M::Mul, k) - A,Bs = first(M.factors), tail(M.factors) + A,Bs = first(M.args), tail(M.args) B = Mul(Bs) ret = zero(eltype(M)) for j = rowsupport(A, k) @@ -177,7 +155,7 @@ function getindex(M::Mul, k) end function getindex(M::Mul, k, j) - A,Bs = first(M.factors), tail(M.factors) + A,Bs = first(M.args), tail(M.args) B = Mul(Bs) ret = zero(eltype(M)) @inbounds for ℓ in (rowsupport(A,k) ∩ colsupport(B,j)) @@ -218,16 +196,16 @@ IndexStyle(::MulArray{<:Any,1}) = IndexLinear() *(A::MulArray, B::Mul) = A.mul * B *(A::Mul, B::MulArray) = A * B.mul -adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.mul.factors))...) -transpose(A::MulArray) = MulArray(reverse(transpose.(A.mul.factors))...) +adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.mul.args))...) +transpose(A::MulArray) = MulArray(reverse(transpose.(A.mul.args))...) struct MulLayout{LAY} <: MemoryLayout layouts::LAY end -MemoryLayout(M::MulArray) = MulLayout(MemoryLayout.(M.mul.factors)) +MemoryLayout(M::MulArray) = MulLayout(MemoryLayout.(M.mul.args)) -_flatten(A::MulArray, B...) = _flatten(A.mul.factors..., B...) -flatten(A::MulArray) = MulArray(Mul(_flatten(A.mul.factors...))) +_flatten(A::MulArray, B...) = _flatten(A.mul.args..., B...) +flatten(A::MulArray) = MulArray(Mul(_flatten(A.mul.args...))) diff --git a/test/applytests.jl b/test/applytests.jl new file mode 100644 index 00000000..d56a0b40 --- /dev/null +++ b/test/applytests.jl @@ -0,0 +1,12 @@ +using LazyArrays, Test +import LazyArrays: materialize, broadcasted, DefaultApplyStyle, Applied + +@test applied(exp,1) isa Applied{DefaultApplyStyle} + +materialize(applied(randn)) + +@test materialize(applied(*, 1)) == 1 + +@test materialize(applied(exp, 1)) === exp(1) +@test materialize(applied(exp, broadcasted(+, 1, 2))) === + materialize(applied(exp, applied(+, 1, 2))) === exp(3) diff --git a/test/multests.jl b/test/multests.jl index 5ac442b6..028ce428 100644 --- a/test/multests.jl +++ b/test/multests.jl @@ -650,9 +650,6 @@ import Base.Broadcast: materialize, materialize! B = materialize(Mul(A,A,A)) @test B isa Matrix{Float64} @test all(B .=== (A*A)*A) - - @test Mul(A,A) * A ≈ A * Mul(A,A) ≈ Mul(A) * Mul(A,A) ≈ A^3 - @test Mul(A,A) * Mul(A,A) ≈ Mul(A) * Mul(A,A,A) ≈ A^4 end @testset "Diagonal and SymTridiagonal" begin diff --git a/test/runtests.jl b/test/runtests.jl index 1e7b3b5d..901bd982 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using Test, LinearAlgebra, LazyArrays, StaticArrays, FillArrays include("memorylayouttests.jl") +include("applytests.jl") include("multests.jl") include("ldivtests.jl") From ac0ae9c25a2d36ef9ebe556faf5cc8226dfa622d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 20 Jan 2019 09:20:17 +0000 Subject: [PATCH 02/15] Ldiv, PInv, and Inv are now Applied --- src/lazyapplying.jl | 54 ++++++++++++++ src/linalg/inv.jl | 158 +++++++++++++++++++---------------------- src/linalg/lazymul.jl | 18 ++++- src/linalg/mul.jl | 7 +- test/applytests.jl | 14 ++-- test/benchmarktests.jl | 6 +- test/ldivtests.jl | 17 ++--- 7 files changed, 167 insertions(+), 107 deletions(-) diff --git a/src/lazyapplying.jl b/src/lazyapplying.jl index 4a34a331..9e95c272 100644 --- a/src/lazyapplying.jl +++ b/src/lazyapplying.jl @@ -17,4 +17,58 @@ struct Applied{Style<:ApplyStyle, F, Args<:Tuple} end applied(f, args...) = Applied(ApplyStyle(f, args...), f, args) +materialize(A::Applied{DefaultApplyStyle,<:Any,<:Tuple{<:Any}}) = + A.f(materialize(first(A.args))) materialize(A::Applied{DefaultApplyStyle}) = A.f(materialize.(A.args)...) + +_apply_eltype(f, a...) = Base.promote_op(f, eltype.(a)...) +eltype(M::Applied) = _apply_eltype(M.f, M.args...) + +similar(M::Applied) = similar(M, eltype(M)) + + +struct ApplyArray{T, N, App<:Applied} <: AbstractArray{T,N} + applied::App +end + +const ApplyVector{T, App<:Applied} = ApplyArray{T, 1, App} +const ApplyMatrix{T, App<:Applied} = ApplyArray{T, 2, App} + + +ApplyArray{T,N}(M::App) where {T,N,App<:Applied} = ApplyArray{T,N,App}(M) +ApplyArray{T}(M::Applied) where {T} = ApplyArray{T,ndims(M)}(M) +ApplyArray(M::Applied) = ApplyArray{eltype(M)}(M) +ApplyVector(M::Applied) = ApplyVector{eltype(M)}(M) +ApplyMatrix(M::Applied) = ApplyMatrix{eltype(M)}(M) + +ApplyArray(f, factors...) = ApplyArray(applied(f, factors...)) +ApplyArray{T}(f, factors...) where T = ApplyArray{T}(applied(f, factors...)) +ApplyArray{T,N}(f, factors...) where {T,N} = ApplyArray{T,N}(applied(f, factors...)) +ApplyVector(f, factors...) = ApplyVector(applied(f, factors...)) +ApplyMatrix(f, factors...) = ApplyMatrix(applied(f, factors...)) + +axes(A::ApplyArray) = axes(A.applied) +size(A::ApplyArray) = map(length, axes(A)) + +IndexStyle(::ApplyArray{<:Any,1}) = IndexLinear() + +@propagate_inbounds getindex(A::ApplyArray{T,N}, kj::Vararg{Int,N}) where {T,N} = + materialize(A.applied)[kj...] + +@inline copyto!(dest::AbstractArray, M::Applied) = _copyto!(MemoryLayout(dest), dest, M) + + +# adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.mul.args))...) +# transpose(A::MulArray) = MulArray(reverse(transpose.(A.mul.args))...) + + +struct ApplyLayout{F, LAY} <: MemoryLayout + f::F + layouts::LAY +end + +MemoryLayout(M::ApplyArray) = ApplyLayout(M.applied.f, MemoryLayout.(M.applied.args)) + + +# _flatten(A::ApplyArray, B...) = _flatten(A.mul.args..., B...) +# flatten(A::MulArray) = MulArray(Mul(_flatten(A.mul.args...))) diff --git a/src/linalg/inv.jl b/src/linalg/inv.jl index 7d2d110f..c90e493e 100644 --- a/src/linalg/inv.jl +++ b/src/linalg/inv.jl @@ -1,66 +1,60 @@ -abstract type AbstractPInv{Style, Typ} end -eltype(::AbstractPInv{<:Any,Typ}) where Typ = eltype(Typ) -eltype(::Type{<:AbstractPInv{<:Any,Typ}}) where Typ = eltype(Typ) -struct PInv{Style, Typ} <: AbstractPInv{Style, Typ} - style::Style - A::Typ -end -struct Inv{Style, Typ} <: AbstractPInv{Style, Typ} - style::Style - A::Typ - function Inv{Style,Typ}(style::Style, A::Typ) where {Style,Typ} - checksquare(A) - new{Style,Typ}(style,A) - end -end +const PInv{Style, Typ} = Applied{<:LayoutApplyStyle{Style}, typeof(pinv), Tuple{Typ}} +const Inv{Style, Typ} = Applied{<:LayoutApplyStyle{Style}, typeof(inv), Tuple{Typ}} -Inv(style::Style, A::Typ) where {Style,Typ} = Inv{Style,Typ}(style, A) +Inv(A) = applied(inv, A) +PInv(A) = applied(pinv, A) -PInv(A) = PInv(MemoryLayout(A), A) -Inv(A) = Inv(MemoryLayout(A), A) +ApplyStyle(::typeof(inv), A::AbstractMatrix) = LayoutApplyStyle((MemoryLayout(A),)) +ApplyStyle(::typeof(pinv), A::AbstractMatrix) = LayoutApplyStyle((MemoryLayout(A),)) +const InvOrPInv = Union{Inv, PInv} +parent(A::PInv) = first(A.args) +parent(A::Inv) = first(A.args) -pinv(A::PInv) = A.A +pinv(A::PInv) = parent(A) function inv(A::PInv) - checksquare(A.A) - A.A + checksquare(parent(A)) + parent(A) end -inv(A::Inv) = A.A +inv(A::Inv) = parent(A) pinv(A::Inv) = inv(A) -parent(A::AbstractPInv) = A.A -size(A::AbstractPInv) = reverse(size(parent(A))) -axes(A::AbstractPInv) = reverse(axes(parent(A))) -size(A::AbstractPInv, k) = size(A)[k] -axes(A::AbstractPInv, k) = axes(A)[k] +size(A::InvOrPInv) = reverse(size(parent(A))) +axes(A::InvOrPInv) = reverse(axes(parent(A))) +size(A::InvOrPInv, k) = size(A)[k] +axes(A::InvOrPInv, k) = axes(A)[k] -abstract type AbstractPInvLayout{ML} <: MemoryLayout end -# struct InverseLayout{ML} <: AbstractPInvLayout{ML} -# layout::ML -# end +const Ldiv{StyleA, StyleB, AType, BType} = + Applied{<:LayoutApplyStyle{<:Tuple{StyleA, StyleB}}, typeof(\), <:Tuple{AType, BType}} -struct PInvLayout{ML} <: AbstractPInvLayout{ML} - layout::ML -end +Ldiv(A, B) = applied(\, A, B) -MemoryLayout(Ai::AbstractPInv) = PInvLayout(MemoryLayout(Ai.A)) +ApplyStyle(::typeof(\), A::AbstractArray, B::AbstractArray) = + LayoutApplyStyle((MemoryLayout(A), MemoryLayout(B))) +size(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractMatrix}) = + (size(L.args[1], 2),size(L.args[2],2)) +size(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) = + (size(L.args[1], 2),) +length(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) = + size(L.args[1], 2) + +struct ArrayLdivArrayStyle{StyleA, StyleB, p, q} <: BroadcastStyle end + +@inline copyto!(dest::AbstractArray, bc::Broadcasted{<:ArrayLdivArrayStyle}) = + _copyto!(MemoryLayout(dest), dest, bc) -const Ldiv{StyleA, StyleB, AType, BType} = - Mul2{<:AbstractPInvLayout{StyleA}, StyleB, <:AbstractPInv{StyleA,AType}, BType} const ArrayLdivArray{styleA, styleB, p, q, T, V} = Ldiv{styleA, styleB, <:AbstractArray{T,p}, <:AbstractArray{V,q}} -const ArrayLdivArrayStyle{StyleA,StyleB,p,q} = - ArrayMulArrayStyle{PInvLayout{StyleA}, StyleB, p, q} const BArrayLdivArray{styleA, styleB, p, q, T, V} = Broadcasted{ArrayLdivArrayStyle{styleA,styleB,p,q}, <:Any, typeof(identity), <:Tuple{<:ArrayLdivArray{styleA,styleB,p,q,T,V}}} @@ -70,24 +64,9 @@ BroadcastStyle(::Type{<:ArrayLdivArray{StyleA,StyleB,p,q}}) where {StyleA,StyleB ArrayLdivArrayStyle{StyleA,StyleB,p,q}() broadcastable(M::ArrayLdivArray) = M -Ldiv(A, B) = Mul(PInv(A), B) -macro lazyldiv(Typ) - esc(quote - LinearAlgebra.ldiv!(A::$Typ, x::AbstractVector) = (x .= LazyArrays.Ldiv(A,x)) - LinearAlgebra.ldiv!(A::$Typ, x::AbstractMatrix) = (x .= LazyArrays.Ldiv(A,x)) - LinearAlgebra.ldiv!(A::$Typ, x::StridedVector) = (x .= LazyArrays.Ldiv(A,x)) - LinearAlgebra.ldiv!(A::$Typ, x::StridedMatrix) = (x .= LazyArrays.Ldiv(A,x)) - - Base.:\(A::$Typ, x::AbstractVector) = PInv(A) * x - Base.:\(A::$Typ, x::AbstractMatrix) = PInv(A) * x - end) -end - -*(A::AbstractPInv, B, C...) = materialize(Mul(A,B, C...)) -*(A::AbstractPInv, B::Mul) = materialize(Mul(A, B.args...)) - -similar(A::AbstractPInv, ::Type{T}) where T = Array{T}(undef, size(A)) +similar(A::InvOrPInv, ::Type{T}) where T = Array{T}(undef, size(A)) +similar(A::Ldiv, ::Type{T}) where T = Array{T}(undef, size(A)) similar(M::ArrayLdivArray, ::Type{T}) where T = Array{T}(undef, size(M)) materialize(M::ArrayLdivArray) = copyto!(similar(M), M) @@ -99,13 +78,13 @@ end if VERSION ≥ v"1.1-pre" function _copyto!(_, dest::AbstractArray, M::ArrayLdivArray) - Ai, B = M.args - ldiv!(dest, factorize(pinv(Ai)), B) + A, B = M.args + ldiv!(dest, factorize(A), B) end else function _copyto!(_, dest::AbstractArray, M::ArrayLdivArray) - Ai, B = M.args - ldiv!(dest, factorize(pinv(Ai)), copy(B)) + A, B = M.args + ldiv!(dest, factorize(A), copy(B)) end end @@ -120,78 +99,87 @@ broadcastable(M::MatLdivVec) = M ### function _copyto!(_, dest::AbstractArray, M::ArrayLdivArray{<:TriangularLayout}) - Ai, B = M.args + A, B = M.args dest ≡ B || (dest .= B) - ldiv!(pinv(Ai), dest) + ldiv!(A, dest) end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{<:TriangularLayout{UPLO,UNIT,<:AbstractColumnMajor}, <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat} - Ai,B = M.args + A,B = M.args B ≡ dest || copyto!(dest, B) - BLAS.trsv!(UPLO, 'N', UNIT, triangulardata(pinv(Ai)), dest) + BLAS.trsv!(UPLO, 'N', UNIT, triangulardata(A), dest) end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{<:TriangularLayout{'U',UNIT,<:AbstractRowMajor}, <:AbstractStridedLayout, T, T}) where {UNIT,T <: BlasFloat} - Ai,B = M.args + A,B = M.args B ≡ dest || copyto!(dest, B) - BLAS.trsv!('L', 'T', UNIT, transpose(triangulardata(pinv(Ai))), dest) + BLAS.trsv!('L', 'T', UNIT, transpose(triangulardata(A)), dest) end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:AbstractRowMajor}, <:AbstractStridedLayout, T, T}) where {UNIT,T <: BlasFloat} - Ai,B = M.args + A,B = M.args B ≡ dest || copyto!(dest, B) - BLAS.trsv!('U', 'T', UNIT, transpose(triangulardata(pinv(Ai))), dest) + BLAS.trsv!('U', 'T', UNIT, transpose(triangulardata(A)), dest) end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{T, <:TriangularLayout{'U',UNIT,<:ConjLayout{<:AbstractRowMajor}}, <:AbstractStridedLayout, T, T}) where {UNIT,T <: BlasFloat} - Ai,B = M.args + A,B = M.args B ≡ dest || copyto!(dest, B) - BLAS.trsv!('L', 'C', UNIT, triangulardata(pinv(Ai))', dest) + BLAS.trsv!('L', 'C', UNIT, triangulardata(A)', dest) end @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T}, M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:ConjLayout{<:AbstractRowMajor}}, <:AbstractStridedLayout, T, T}) where {UNIT,T <: BlasFloat} - Ai,B = M.args + A,B = M.args B ≡ dest || copyto!(dest, B) - BLAS.trsv!('U', 'C', UNIT, triangulardata(pinv(Ai))', dest) + BLAS.trsv!('U', 'C', UNIT, triangulardata(A)', dest) end function _copyto!(_, dest::AbstractMatrix, M::MatLdivMat{<:TriangularLayout}) A,X = M.args size(dest,2) == size(X,2) || thow(DimensionMismatch("Dimensions must match")) @views for j in axes(dest,2) - dest[:,j] .= Mul(A, X[:,j]) + dest[:,j] .= Ldiv(A, X[:,j]) end dest end -function _PInvMatrix end +const PInvMatrix{T,App<:PInv} = ApplyMatrix{T,App} +const InvMatrix{T,App<:Inv} = ApplyMatrix{T,App} -struct PInvMatrix{T, PINV<:AbstractPInv} <: AbstractMatrix{T} - pinv::PINV - global _PInvMatrix(pinv::PINV) where PINV = new{eltype(pinv), PINV}(pinv) +PInvMatrix(A) = ApplyMatrix(pinv, A) +function InvMatrix(A) + checksquare(A) + ApplyMatrix(inv, A) end -const InvMatrix{T, INV<:Inv} = PInvMatrix{T,INV} - -PInvMatrix(M) = _PInvMatrix(PInv(M)) -InvMatrix(M) = _PInvMatrix(Inv(M)) - -axes(A::PInvMatrix) = axes(A.pinv) +axes(A::PInvMatrix) = reverse(axes(parent(A.applied))) size(A::PInvMatrix) = map(length, axes(A)) @propagate_inbounds getindex(A::PInvMatrix{T}, k::Int, j::Int) where T = - (A.pinv*[Zeros(j-1); one(T); Zeros(size(A,2) - j)])[k] + (parent(A.applied)\[Zeros(j-1); one(T); Zeros(size(A,2) - j)])[k] + +@propagate_inbounds getindex(A::InvMatrix{T}, k::Int, j::Int) where T = + (parent(A.applied)\[Zeros(j-1); one(T); Zeros(size(A,2) - j)])[k] -MemoryLayout(M::PInvMatrix) = MemoryLayout(M.pinv) + +@inline function _copyto!(_, dest::AbstractArray, M::MatMulVec{<:ApplyLayout{typeof(inv)}}) + Ai,b = M.args + dest .= Ldiv(parent(Ai.applied), b) +end + +@inline function _copyto!(_, dest::AbstractArray, M::MatMulVec{<:ApplyLayout{typeof(pinv)}}) + Ai,b = M.args + dest .= Ldiv(parent(Ai.applied), b) +end diff --git a/src/linalg/lazymul.jl b/src/linalg/lazymul.jl index d055eea5..6119a5ac 100644 --- a/src/linalg/lazymul.jl +++ b/src/linalg/lazymul.jl @@ -65,5 +65,21 @@ macro lazylmul(Typ) end) end +macro lazyldiv(Typ) + esc(quote + LinearAlgebra.ldiv!(A::$Typ, x::AbstractVector) = (x .= LazyArrays.Ldiv(A,x)) + LinearAlgebra.ldiv!(A::$Typ, x::AbstractMatrix) = (x .= LazyArrays.Ldiv(A,x)) + LinearAlgebra.ldiv!(A::$Typ, x::StridedVector) = (x .= LazyArrays.Ldiv(A,x)) + LinearAlgebra.ldiv!(A::$Typ, x::StridedMatrix) = (x .= LazyArrays.Ldiv(A,x)) + + Base.:\(A::$Typ, x::AbstractVector) = LazyArrays.materialize(LazyArrays.Ldiv(A,x)) + Base.:\(A::$Typ, x::AbstractMatrix) = LazyArrays.materialize(LazyArrays.Ldiv(A,x)) + end) +end -@lazymul AddMatrix +@lazymul ApplyArray +@lazylmul ApplyArray +@lazyldiv ApplyArray +@lazymul BroadcastArray +@lazylmul BroadcastArray +@lazyldiv BroadcastArray diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index d89433ef..c3af04a3 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -15,9 +15,6 @@ Mul(A...) = applied(*, A...) const Mul2{StyleA, StyleB, AType, BType} = Mul{<:Tuple{StyleA,StyleB}, <:Tuple{AType,BType}} -_mul_eltype(a) = eltype(a) -_mul_eltype(a, b...) = Base.promote_op(*, eltype(a), _mul_eltype(b...)) -eltype(M::Mul) = _mul_eltype(M.args...) size(M::Mul, p::Int) = size(M)[p] axes(M::Mul, p::Int) = axes(M)[p] ndims(M::Mul) = ndims(last(M.args)) @@ -31,7 +28,7 @@ _mul_axes(ax1, (_,ax2)::Tuple{<:Any,<:Any}) = (ax1,ax2) axes(M::Mul) = _mul_axes(axes(first(M.args),1), axes(last(M.args))) axes(M::Mul{Tuple{}}) = () -similar(M::Mul) = similar(M, eltype(M)) + @@ -80,8 +77,6 @@ _materialize(M::Mul, _) = lmaterialize(M) _materialize(M::Mul2, _) = error("Cannot materialize $M") materialize(M::Mul) = _materialize(M, axes(M)) -@inline copyto!(dest::AbstractArray, M::Mul) = _copyto!(MemoryLayout(dest), dest, M) - #### # Matrix * Vector diff --git a/test/applytests.jl b/test/applytests.jl index d56a0b40..1e0ab058 100644 --- a/test/applytests.jl +++ b/test/applytests.jl @@ -1,12 +1,14 @@ using LazyArrays, Test import LazyArrays: materialize, broadcasted, DefaultApplyStyle, Applied -@test applied(exp,1) isa Applied{DefaultApplyStyle} +@testset "Applied" begin + @test applied(exp,1) isa Applied{DefaultApplyStyle} -materialize(applied(randn)) + @test materialize(applied(randn)) isa Float64 -@test materialize(applied(*, 1)) == 1 + @test materialize(applied(*, 1)) == 1 -@test materialize(applied(exp, 1)) === exp(1) -@test materialize(applied(exp, broadcasted(+, 1, 2))) === - materialize(applied(exp, applied(+, 1, 2))) === exp(3) + @test materialize(applied(exp, 1)) === exp(1) + @test materialize(applied(exp, broadcasted(+, 1, 2))) === + materialize(applied(exp, applied(+, 1, 2))) === exp(3) +end diff --git a/test/benchmarktests.jl b/test/benchmarktests.jl index 691c984d..4d638941 100644 --- a/test/benchmarktests.jl +++ b/test/benchmarktests.jl @@ -3,6 +3,10 @@ using LazyArrays, BenchmarkTools +@testset "Applied" begin + @test @belapsed(materialize(applied(exp, $x))) ≤ 2(@belapsed exp($x)) +end + @testset "concat" begin A = Vcat(Vector(1:10), Vector(1:20)) b = Array{Int}(undef, 30) @@ -19,4 +23,4 @@ using LazyArrays, BenchmarkTools A = Hcat(1:10, 2:11) b = Array{Int}(undef, 10, 2) @test @belapsed(copyto!($b,$A)) < @belapsed(hcat($A.arrays...)) -end +en diff --git a/test/ldivtests.jl b/test/ldivtests.jl index 992fdb3e..8e508623 100644 --- a/test/ldivtests.jl +++ b/test/ldivtests.jl @@ -1,12 +1,12 @@ using LazyArrays, LinearAlgebra, Test -import LazyArrays: ArrayLdivArrayStyle +import LazyArrays: ArrayLdivArrayStyle, InvMatrix import Base.Broadcast: materialize @testset "Ldiv" begin @testset "Float64 \\ *" begin A = randn(5,5) b = randn(5) - M = Mul(Inv(A),b) + M = Ldiv(A,b) @test size(M) == (5,) @test similar(M) isa Vector{Float64} @@ -15,14 +15,15 @@ import Base.Broadcast: materialize @test Base.BroadcastStyle(typeof(Ldiv(A,b))) isa ArrayLdivArrayStyle @test all(copyto!(similar(b), Ldiv(A,b)) .=== - (similar(b) .= Ldiv(A,b)) .=== Inv(A) * b .=== + (similar(b) .= Ldiv(A,b)) .=== InvMatrix(A) * b .=== materialize(Ldiv(A,b)) .=== (A\b) .=== (b̃ = copy(b); LAPACK.gesv!(copy(A), b̃); b̃)) @test copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) ≈ UpperTriangular(A) \ b @test all(copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) .=== - (similar(b) .= Ldiv(UpperTriangular(A),b)) .=== Inv(UpperTriangular(A))*b .=== + (similar(b) .= Ldiv(UpperTriangular(A),b)) .=== + InvMatrix(UpperTriangular(A))*b .=== BLAS.trsv('U', 'N', 'N', A, b) ) @@ -34,7 +35,7 @@ import Base.Broadcast: materialize BLAS.trsv('U', 'T', 'N', A, b)) b = randn(5) + im*randn(5) - @test Inv(A) * b ≈ Matrix(A) \ b + @test InvMatrix(A) * b ≈ Matrix(A) \ b end @@ -63,7 +64,7 @@ import Base.Broadcast: materialize BLAS.trsv('U', 'T', 'N', A, b)) b = randn(5) - @test Inv(A) * b ≈ A \ b + @test InvMatrix(A) * b ≈ A \ b end @testset "Triangular \\ matrix" begin @@ -80,8 +81,8 @@ import Base.Broadcast: materialize A_orig = copy(A) b = randn(5) b_orig = copy(b) - @test_throws DimensionMismatch Inv(A) - @test PInv(A) * b == (A\b) + @test_throws DimensionMismatch InvMatrix(A) + @test PInvMatrix(A) * b == (A\b) @test all(b .=== b_orig) end From b8217b7eb9869ef297fa726da9bb17cece534ccc Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 20 Jan 2019 16:45:33 +0000 Subject: [PATCH 03/15] Add is now Applied(+, ...) --- src/lazyapplying.jl | 32 ++++++++++++++++++++++++++++++-- src/lazybroadcasting.jl | 11 ----------- src/linalg/blasmul.jl | 10 ++++++++-- src/linalg/inv.jl | 2 ++ src/linalg/mul.jl | 1 + test/memorylayouttests.jl | 11 +++++++++-- test/multests.jl | 13 +++++++------ test/runtests.jl | 1 + test/setoptests.jl | 6 +++--- 9 files changed, 61 insertions(+), 26 deletions(-) diff --git a/src/lazyapplying.jl b/src/lazyapplying.jl index 9e95c272..c86a90a0 100644 --- a/src/lazyapplying.jl +++ b/src/lazyapplying.jl @@ -21,8 +21,7 @@ materialize(A::Applied{DefaultApplyStyle,<:Any,<:Tuple{<:Any}}) = A.f(materialize(first(A.args))) materialize(A::Applied{DefaultApplyStyle}) = A.f(materialize.(A.args)...) -_apply_eltype(f, a...) = Base.promote_op(f, eltype.(a)...) -eltype(M::Applied) = _apply_eltype(M.f, M.args...) + similar(M::Applied) = similar(M, eltype(M)) @@ -44,6 +43,7 @@ ApplyMatrix(M::Applied) = ApplyMatrix{eltype(M)}(M) ApplyArray(f, factors...) = ApplyArray(applied(f, factors...)) ApplyArray{T}(f, factors...) where T = ApplyArray{T}(applied(f, factors...)) ApplyArray{T,N}(f, factors...) where {T,N} = ApplyArray{T,N}(applied(f, factors...)) + ApplyVector(f, factors...) = ApplyVector(applied(f, factors...)) ApplyMatrix(f, factors...) = ApplyMatrix(applied(f, factors...)) @@ -72,3 +72,31 @@ MemoryLayout(M::ApplyArray) = ApplyLayout(M.applied.f, MemoryLayout.(M.applied.a # _flatten(A::ApplyArray, B...) = _flatten(A.mul.args..., B...) # flatten(A::MulArray) = MulArray(Mul(_flatten(A.mul.args...))) + + + +const Add{Factors<:Tuple} = Applied{<:Any, typeof(+), Factors} + +size(M::Add, p::Int) = size(M)[p] +axes(M::Add, p::Int) = axes(M)[p] +ndims(M::Add) = ndims(first(M.args)) + +length(M::Add) = prod(size(M)) +size(M::Add) = length.(axes(M)) +axes(M::Add) = axes(first(M.args)) + + +eltype(M::Add) = Base._return_type(+, eltype.(M.args)) + +const AddArray{T,N,Factors<:Tuple} = ApplyArray{T,N,<:Add{Factors}} +const AddVector{T,Factors<:Tuple} = AddArray{T,1,Factors} +const AddMatrix{T,Factors<:Tuple} = AddArray{T,2,Factors} + +AddArray(factors...) = ApplyArray(+, factors...) + +""" + Add(A1, A2, …, AN) + +A lazy representation of `A1 + A2 + … + AN`; i.e., a shorthand for `applied(+, A1, A2, …, AN)`. +""" +Add(As...) = applied(+, As...) diff --git a/src/lazybroadcasting.jl b/src/lazybroadcasting.jl index 8d94ed5b..a2387d13 100644 --- a/src/lazybroadcasting.jl +++ b/src/lazybroadcasting.jl @@ -90,14 +90,3 @@ broadcasted(::LazyArrayStyle{N}, ::typeof(*), a::AbstractArray{T,N}, b::Zeros{V, broadcast(DefaultArrayStyle{N}(), *, a, b) broadcasted(::LazyArrayStyle{N}, ::typeof(*), a::Zeros{T,N}, b::AbstractArray{V,N}) where {T,V,N} = broadcast(DefaultArrayStyle{N}(), *, a, b) - -const Add = BroadcastArray{<:Any, <:Any, <:Broadcasted{<:Any, <:Any, typeof(+)}} -const AddVector = Add{<:Any, 1} -const AddMatrix = Add{<:Any, 2} - -""" - Add(A1, A2, …, AN) - -A lazy representation of `A1 .+ A2 .+ … .+ AN`; i.e., a shorthand for `BroadcastArray(+, A1, A2, …, AN)`. -""" -Add(As...) = BroadcastArray(+, As...) diff --git a/src/linalg/blasmul.jl b/src/linalg/blasmul.jl index 79908d9f..d1abf8f5 100644 --- a/src/linalg/blasmul.jl +++ b/src/linalg/blasmul.jl @@ -184,6 +184,9 @@ end function materialize!(M::MatMulMatAdd) α, A, B, β, C = M.α, M.A, M.B, M.β, M.C + if C ≡ B + B = copy(B) + end ts = tile_size(eltype(A), eltype(B), eltype(C)) if iszero(β) # false is a "strong" zero to wipe out NaNs ts == 0 ? default_blasmul!(α, A, B, false, C) : tiled_blasmul!(ts, α, A, B, false, C) @@ -194,19 +197,22 @@ end function materialize!(M::MatMulVecAdd) α, A, B, β, C = M.α, M.A, M.B, M.β, M.C + if C ≡ B + B = copy(B) + end default_blasmul!(α, A, B, iszero(β) ? false : β, C) end for MulAdd_ in [MatMulMatAdd, MatMulVecAdd] # `MulAdd{<:BroadcastLayout{typeof(+)}}` cannot "win" against # `MatMulMatAdd` and `MatMulVecAdd` hence `@eval`: - @eval function materialize!(M::$MulAdd_{<:BroadcastLayout{typeof(+)}}) + @eval function materialize!(M::$MulAdd_{<:ApplyLayout{typeof(+)}}) α, A, B, β, C = M.α, M.A, M.B, M.β, M.C if C ≡ B B = copy(B) end lmul!(β, C) - for A in A.broadcasted.args + for A in A.applied.args C .= α .* Mul(A, B) .+ C end C diff --git a/src/linalg/inv.jl b/src/linalg/inv.jl index c90e493e..6e7ec9f0 100644 --- a/src/linalg/inv.jl +++ b/src/linalg/inv.jl @@ -30,6 +30,7 @@ size(A::InvOrPInv) = reverse(size(parent(A))) axes(A::InvOrPInv) = reverse(axes(parent(A))) size(A::InvOrPInv, k) = size(A)[k] axes(A::InvOrPInv, k) = axes(A)[k] +eltype(A::InvOrPInv) = eltype(parent(A)) @@ -47,6 +48,7 @@ size(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) = (size(L.args[1], 2),) length(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) = size(L.args[1], 2) +eltype(M::Ldiv) = promote_type(eltype.(M.args)...) struct ArrayLdivArrayStyle{StyleA, StyleB, p, q} <: BroadcastStyle end diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index c3af04a3..03e9c049 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -21,6 +21,7 @@ ndims(M::Mul) = ndims(last(M.args)) length(M::Mul) = prod(size(M)) size(M::Mul) = length.(axes(M)) +eltype(M::Mul) = Base.promote_op(*, eltype.(M.args)...) _mul_axes(ax1, ::Tuple{}) = (ax1,) _mul_axes(ax1, ::Tuple{<:Any}) = (ax1,) diff --git a/test/memorylayouttests.jl b/test/memorylayouttests.jl index 1b9c5224..55d1f9fe 100644 --- a/test/memorylayouttests.jl +++ b/test/memorylayouttests.jl @@ -5,7 +5,7 @@ using LazyArrays, LinearAlgebra, FillArrays, Test UnitUpperTriangularLayout, LowerTriangularLayout, UnitLowerTriangularLayout, ScalarLayout, hermitiandata, symmetricdata, FillLayout, ZerosLayout, - VcatLayout, BroadcastLayout, Add + VcatLayout, BroadcastLayout, Add, AddArray, ApplyLayout struct FooBar end struct FooNumber <: Number end @@ -145,7 +145,14 @@ struct FooNumber <: Number end @testset "BroadcastArray" begin A = [1.0 2; 3 4] - @test MemoryLayout(Add(A, Fill(0, (2, 2)), Zeros(2, 2))) == + @test MemoryLayout(BroadcastArray(+, A, Fill(0, (2, 2)), Zeros(2, 2))) == BroadcastLayout(+, (DenseColumnMajor(), FillLayout(), ZerosLayout())) end + + @testset "ApplyArray" begin + A = [1.0 2; 3 4] + @test eltype(AddArray(A, Fill(0, (2, 2)), Zeros(2, 2))) == Float64 + @test MemoryLayout(AddArray(A, Fill(0, (2, 2)), Zeros(2, 2))) == + ApplyLayout(+, (DenseColumnMajor(), FillLayout(), ZerosLayout())) + end end diff --git a/test/multests.jl b/test/multests.jl index 028ce428..85230262 100644 --- a/test/multests.jl +++ b/test/multests.jl @@ -1,7 +1,8 @@ using Test, LinearAlgebra, LazyArrays, StaticArrays, FillArrays -import LazyArrays: MulAdd, MemoryLayout, DenseColumnMajor, DiagonalLayout, SymTridiagonalLayout, Add +import LazyArrays: MulAdd, MemoryLayout, DenseColumnMajor, DiagonalLayout, SymTridiagonalLayout, Add, AddArray import Base.Broadcast: materialize, materialize! + @testset "Mul" begin @testset "eltype" begin @inferred(eltype(Mul(zeros(Int,2,2), zeros(Float64,2)))) == Float64 @@ -640,7 +641,7 @@ import Base.Broadcast: materialize, materialize! Ac = A' blasnoalloc(c, 2.0, Ac, x, 3.0, y) @test @allocated(blasnoalloc(c, 2.0, Ac, x, 3.0, y)) == 0 - Aa = Add(A, Ac) + Aa = AddArray(A, Ac) blasnoalloc(c, 2.0, Aa, x, 3.0, y) @test_broken @allocated(blasnoalloc(c, 2.0, Aa, x, 3.0, y)) == 0 end @@ -674,8 +675,8 @@ end @testset "Add" begin @testset "gemv Float64" begin - for A in (Add(randn(5,5), randn(5,5)), - Add(randn(5,5), view(randn(9, 5), 1:2:9, :))), + for A in (AddArray(randn(5,5), randn(5,5)), + AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))), b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)) à = copy(A) @@ -718,8 +719,8 @@ end end @testset "gemm" begin - for A in (Add(randn(5,5), randn(5,5)), - Add(randn(5,5), view(randn(9, 5), 1:2:9, :))), + for A in (AddArray(randn(5,5), randn(5,5)), + AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))), B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)) diff --git a/test/runtests.jl b/test/runtests.jl index 901bd982..ac9ccdef 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ include("memorylayouttests.jl") include("applytests.jl") include("multests.jl") include("ldivtests.jl") +include("setoptests.jl") @testset "concat" begin diff --git a/test/setoptests.jl b/test/setoptests.jl index 5502c256..1acba670 100644 --- a/test/setoptests.jl +++ b/test/setoptests.jl @@ -1,6 +1,6 @@ using LazyArrays, Test using LazyArrays.SetOperations -import LazyArrays.SetOperations: ∪, ∩, unioned, intersected, setdiffed +import LazyArrays.SetOperations: unioned, intersected, setdiffed @testset "Set operations" begin @test SetOperations.Style(Int) == SetOperations.VectorStyle() @@ -10,8 +10,8 @@ import LazyArrays.SetOperations: ∪, ∩, unioned, intersected, setdiffed @test 3 ∉ unioned(1,2.0) @test issubset([1,2], unioned(1,2.0)) - @test 1 ∪ 2.0 == [1.0,2.0] - @test [1,2] ∪ Set([2,4]) == Set([2,4]) ∪ [1,2] == SetOperations.union(Set([2,4]) , [1,2]) == + @test SetOperations.:∪(1, 2.0) == [1.0,2.0] + @test SetOperations.:∪([1,2], Set([2,4])) == Set([2,4]) ∪ [1,2] == SetOperations.union(Set([2,4]) , [1,2]) == SetOperations.union([1,2], Set([2,4])) == Set([1,2,4]) @test Base.union([1,2], Set([2,4])) == Base.:∪([1,2], Set([2,4])) == [1,2,4] From 10feb38a12137df55277e6b9e77921cb449f62e9 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 20 Jan 2019 16:51:26 +0000 Subject: [PATCH 04/15] MulArray is now an ApplyArray --- src/linalg/mul.jl | 24 +++++------------------- 1 file changed, 5 insertions(+), 19 deletions(-) diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index 03e9c049..86d37b23 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -161,9 +161,7 @@ function getindex(M::Mul, k, j) end -struct MulArray{T, N, MUL<:Mul} <: AbstractArray{T,N} - mul::MUL -end +const MulArray{T, N, MUL<:Mul} = ApplyArray{T, N, MUL} const MulVector{T, MUL<:Mul} = MulArray{T, 1, MUL} const MulMatrix{T, MUL<:Mul} = MulArray{T, 2, MUL} @@ -181,26 +179,14 @@ MulArray{T,N}(factors...) where {T,N} = MulArray{T,N}(Mul(factors...)) MulVector(factors...) = MulVector(Mul(factors...)) MulMatrix(factors...) = MulMatrix(Mul(factors...)) -axes(A::MulArray) = axes(A.mul) -size(A::MulArray) = map(length, axes(A)) - IndexStyle(::MulArray{<:Any,1}) = IndexLinear() -@propagate_inbounds getindex(A::MulArray, kj::Int...) = A.mul[kj...] - -*(A::MulArray, B::MulArray) = A.mul * B.mul -*(A::MulArray, B::Mul) = A.mul * B -*(A::Mul, B::MulArray) = A * B.mul +@propagate_inbounds getindex(A::MulArray, kj::Int...) = A.applied[kj...] -adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.mul.args))...) -transpose(A::MulArray) = MulArray(reverse(transpose.(A.mul.args))...) - - -struct MulLayout{LAY} <: MemoryLayout - layouts::LAY -end +*(A::MulArray, B::MulArray) = MulArray(A, B) -MemoryLayout(M::MulArray) = MulLayout(MemoryLayout.(M.mul.args)) +adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.applied.args))...) +transpose(A::MulArray) = MulArray(reverse(transpose.(A.applied.args))...) _flatten(A::MulArray, B...) = _flatten(A.mul.args..., B...) From 3574e6720bc8cc4e3ae902be4538fb1885a50c33 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 20 Jan 2019 17:38:43 +0000 Subject: [PATCH 05/15] ApplyArray tests --- src/lazyapplying.jl | 17 +++++++++++++++++ src/linalg/mul.jl | 5 ++++- test/applytests.jl | 21 ++++++++++++++++++++- 3 files changed, 41 insertions(+), 2 deletions(-) diff --git a/src/lazyapplying.jl b/src/lazyapplying.jl index c86a90a0..b25593de 100644 --- a/src/lazyapplying.jl +++ b/src/lazyapplying.jl @@ -26,6 +26,23 @@ materialize(A::Applied{DefaultApplyStyle}) = A.f(materialize.(A.args)...) similar(M::Applied) = similar(M, eltype(M)) +struct MatrixFunctionStyle{F} <: ApplyStyle end + +for f in (:exp, :sin, :cos, :sqrt) + @eval ApplyStyle(::typeof($f), ::AbstractMatrix) = MatrixFunctionStyle{typeof($f)}() +end + +materialize(A::Applied{<:MatrixFunctionStyle,<:Any,<:Tuple{<:Any}}) = + A.f(materialize(first(A.args))) + +axes(A::Applied{<:MatrixFunctionStyle}) = axes(first(A.args)) +size(A::Applied{<:MatrixFunctionStyle}) = size(first(A.args)) +eltype(A::Applied{<:MatrixFunctionStyle}) = eltype(first(A.args)) + +getindex(A::Applied{<:MatrixFunctionStyle}, k::Int, j::Int) = + materialize(A)[k,j] + + struct ApplyArray{T, N, App<:Applied} <: AbstractArray{T,N} applied::App end diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index 86d37b23..789b946a 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -150,9 +150,12 @@ function getindex(M::Mul, k) ret end +_mul(A) = A +_mul(A,B,C...) = Mul(A,B,C...) + function getindex(M::Mul, k, j) A,Bs = first(M.args), tail(M.args) - B = Mul(Bs) + B = _mul(Bs...) ret = zero(eltype(M)) @inbounds for ℓ in (rowsupport(A,k) ∩ colsupport(B,j)) ret += A[k,ℓ] * B[ℓ,j] diff --git a/test/applytests.jl b/test/applytests.jl index 1e0ab058..6773f27a 100644 --- a/test/applytests.jl +++ b/test/applytests.jl @@ -1,5 +1,6 @@ using LazyArrays, Test -import LazyArrays: materialize, broadcasted, DefaultApplyStyle, Applied +import LazyArrays: materialize, broadcasted, DefaultApplyStyle, Applied, + ApplyArray, ApplyMatrix, ApplyVector @testset "Applied" begin @test applied(exp,1) isa Applied{DefaultApplyStyle} @@ -12,3 +13,21 @@ import LazyArrays: materialize, broadcasted, DefaultApplyStyle, Applied @test materialize(applied(exp, broadcasted(+, 1, 2))) === materialize(applied(exp, applied(+, 1, 2))) === exp(3) end + +@testset "ApplyArray" begin + A = randn(2,2) + M = ApplyMatrix(exp, A) + @test eltype(M) == Float64 + @test M == exp(A) + + b = randn(2) + c = MulVector(ApplyMatrix(exp, A), b) + + @test axes(c) == (Base.OneTo(2),) + + @test c[1] == c[1,1] + @test exp(A)*b == c + + @test ApplyArray(+,[1,2],[3,4]) == ApplyVector(+,[1,2],[3,4]) == + ApplyArray(+,[1,2],[3,4]) +end From 682b7f5cc9d9a86164df2308e6d873ad9cec88eb Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 21 Jan 2019 21:42:26 +0000 Subject: [PATCH 06/15] add tests --- test/applytests.jl | 2 +- test/multests.jl | 17 +++++++++++++++++ 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/test/applytests.jl b/test/applytests.jl index 6773f27a..72a02292 100644 --- a/test/applytests.jl +++ b/test/applytests.jl @@ -26,7 +26,7 @@ end @test axes(c) == (Base.OneTo(2),) @test c[1] == c[1,1] - @test exp(A)*b == c + @test exp(A)*b ≈ c @test ApplyArray(+,[1,2],[3,4]) == ApplyVector(+,[1,2],[3,4]) == ApplyArray(+,[1,2],[3,4]) diff --git a/test/multests.jl b/test/multests.jl index 85230262..37e04991 100644 --- a/test/multests.jl +++ b/test/multests.jl @@ -670,6 +670,15 @@ import Base.Broadcast: materialize, materialize! M = MulArray(A,A) @test Matrix(M) ≈ A^2 end + + @testset "#14" begin + A = ones(1,1) * 1e200 + B = ones(1,1) * 1e150 + C = ones(1,1) * 1e-300 + + @test materialize(Mul(A, Mul(B,C))) == A*(B*C) + @test materialize(Mul(A , Mul(B , C), C)) == A * (B*C) * C + end end @@ -761,3 +770,11 @@ end end end end + + +N = 2 +A = randn(N,N); B = randn(N,N); C = randn(N,N); R1 = similar(A); R2 = similar(A) + +Mul(A, Mul(B, C)) + +R2 .= Mul(A, Mul(B, C)) From e8492accfabd857013e13d8c4e1e26c7c1ba7fa9 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 22 Jan 2019 07:16:01 +0000 Subject: [PATCH 07/15] Work on #15 --- src/lazyapplying.jl | 2 ++ src/linalg/blasbroadcasting.jl | 2 -- src/linalg/blasmul.jl | 3 +-- src/linalg/inv.jl | 2 -- src/linalg/mul.jl | 13 +++++++++---- test/multests.jl | 17 +++++++++-------- 6 files changed, 21 insertions(+), 18 deletions(-) diff --git a/src/lazyapplying.jl b/src/lazyapplying.jl index b25593de..17b78aa6 100644 --- a/src/lazyapplying.jl +++ b/src/lazyapplying.jl @@ -73,7 +73,9 @@ IndexStyle(::ApplyArray{<:Any,1}) = IndexLinear() materialize(A.applied)[kj...] @inline copyto!(dest::AbstractArray, M::Applied) = _copyto!(MemoryLayout(dest), dest, M) +@inline _copyto!(_, dest::AbstractArray, M::Applied) = copyto!(dest, materialize(M)) +broadcastable(M::Applied) = M # adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.mul.args))...) # transpose(A::MulArray) = MulArray(reverse(transpose.(A.mul.args))...) diff --git a/src/linalg/blasbroadcasting.jl b/src/linalg/blasbroadcasting.jl index 28b46fbc..954da474 100644 --- a/src/linalg/blasbroadcasting.jl +++ b/src/linalg/blasbroadcasting.jl @@ -55,8 +55,6 @@ similar(M::Broadcasted{<:ArrayMulArrayStyle}, ::Type{ElType}) where ElType = Array{ElType}(undef,size(M.args[1])) - -broadcastable(M::ArrayMulArray) = M instantiate(bc::Broadcasted{<:ArrayMulArrayStyle}) = bc diff --git a/src/linalg/blasmul.jl b/src/linalg/blasmul.jl index d1abf8f5..d2c57dad 100644 --- a/src/linalg/blasmul.jl +++ b/src/linalg/blasmul.jl @@ -42,8 +42,7 @@ BroadcastStyle(::Type{<:MatMulVecAdd{StyleA,StyleB,StyleC}}) where {StyleA,Style BroadcastStyle(::Type{<:MatMulMatAdd{StyleA,StyleB,StyleC}}) where {StyleA,StyleB,StyleC} = ArrayMulArrayStyle{StyleA,StyleB,2,2}() -broadcastable(M::MatMulMatAdd) = M -broadcastable(M::MatMulVecAdd) = M +broadcastable(M::MulAdd) = M const BlasMatMulVec{StyleA,StyleB,StyleC,T} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractMatrix{T},<:AbstractVector{T},<:AbstractVector{T}} const BlasMatMulMat{StyleA,StyleB,StyleC,T} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractMatrix{T},<:AbstractMatrix{T},<:AbstractMatrix{T}} diff --git a/src/linalg/inv.jl b/src/linalg/inv.jl index 6e7ec9f0..b6c06a2f 100644 --- a/src/linalg/inv.jl +++ b/src/linalg/inv.jl @@ -64,7 +64,6 @@ const BArrayLdivArray{styleA, styleB, p, q, T, V} = BroadcastStyle(::Type{<:ArrayLdivArray{StyleA,StyleB,p,q}}) where {StyleA,StyleB,p,q} = ArrayLdivArrayStyle{StyleA,StyleB,p,q}() -broadcastable(M::ArrayLdivArray) = M similar(A::InvOrPInv, ::Type{T}) where T = Array{T}(undef, size(A)) @@ -93,7 +92,6 @@ end const MatLdivVec{styleA, styleB, T, V} = ArrayLdivArray{styleA, styleB, 2, 1, T, V} const MatLdivMat{styleA, styleB, T, V} = ArrayLdivArray{styleA, styleB, 2, 2, T, V} -broadcastable(M::MatLdivVec) = M ### diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index 789b946a..4d7bebcc 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -11,7 +11,7 @@ const Mul{Styles<:Tuple, Factors<:Tuple} = Applied{<:LayoutApplyStyle{Styles}, t ApplyStyle(::typeof(*), args::AbstractArray...) = LayoutApplyStyle(MemoryLayout.(args)) -Mul(A...) = applied(*, A...) +Mul(args...) = Applied(LayoutApplyStyle(MemoryLayout.(args)), *, args) const Mul2{StyleA, StyleB, AType, BType} = Mul{<:Tuple{StyleA,StyleB}, <:Tuple{AType,BType}} @@ -19,6 +19,11 @@ size(M::Mul, p::Int) = size(M)[p] axes(M::Mul, p::Int) = axes(M)[p] ndims(M::Mul) = ndims(last(M.args)) +_mul_ndims(::Type{Tuple{A}}) where A = ndims(A) +_mul_ndims(::Type{Tuple{A,B}}) where {A,B} = ndims(B) +ndims(::Type{<:Mul{<:Any,Args}}) where Args = _mul_ndims(Args) + + length(M::Mul) = prod(size(M)) size(M::Mul) = length.(axes(M)) eltype(M::Mul) = Base.promote_op(*, eltype.(M.args)...) @@ -75,7 +80,7 @@ similar(M::ArrayMuls, ::Type{T}) where T = similar(M, T, axes(M)) _materialize(M::ArrayMulArray, _) = copyto!(similar(M), M) _materialize(M::ArrayMuls, _) = lmaterialize(M) _materialize(M::Mul, _) = lmaterialize(M) -_materialize(M::Mul2, _) = error("Cannot materialize $M") +_materialize(M::Mul2, _) = *(materialize.(M.args)...) materialize(M::Mul) = _materialize(M, axes(M)) @@ -140,7 +145,7 @@ getindex(M::MatMulMat, kj::CartesianIndex{2}) = M[kj[1], kj[2]] # MulArray ##### -function getindex(M::Mul, k) +function getindex(M::Mul, k::Integer) A,Bs = first(M.args), tail(M.args) B = Mul(Bs) ret = zero(eltype(M)) @@ -153,7 +158,7 @@ end _mul(A) = A _mul(A,B,C...) = Mul(A,B,C...) -function getindex(M::Mul, k, j) +function getindex(M::Mul, k::Integer, j::Integer) A,Bs = first(M.args), tail(M.args) B = _mul(Bs...) ret = zero(eltype(M)) diff --git a/test/multests.jl b/test/multests.jl index 37e04991..70429f75 100644 --- a/test/multests.jl +++ b/test/multests.jl @@ -679,6 +679,15 @@ import Base.Broadcast: materialize, materialize! @test materialize(Mul(A, Mul(B,C))) == A*(B*C) @test materialize(Mul(A , Mul(B , C), C)) == A * (B*C) * C end + + @testset "#15" begin + N = 2 + A = randn(N,N); B = randn(N,N); C = randn(N,N); R1 = similar(A); R2 = similar(A) + M = Mul(A, Mul(B, C)) + @test ndims(M) == ndims(typeof(M)) == 2 + @test eltype(M) == Float64 + @test_broken all(copyto!(R1, M) .=== A*(B*C) .=== (R2 .= M)) + end end @@ -770,11 +779,3 @@ end end end end - - -N = 2 -A = randn(N,N); B = randn(N,N); C = randn(N,N); R1 = similar(A); R2 = similar(A) - -Mul(A, Mul(B, C)) - -R2 .= Mul(A, Mul(B, C)) From ac0f5d1178ea903c936eaab51e7f531938dd1645 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 22 Jan 2019 21:33:45 +0000 Subject: [PATCH 08/15] Seperate add to own file --- src/lazyapplying.jl | 28 ------------------------- src/linalg/add.jl | 48 +++++++++++++++++++++++++++++++++++++++++++ src/linalg/blasmul.jl | 16 --------------- src/linalg/linalg.jl | 1 + 4 files changed, 49 insertions(+), 44 deletions(-) create mode 100644 src/linalg/add.jl diff --git a/src/lazyapplying.jl b/src/lazyapplying.jl index 17b78aa6..09fd8f01 100644 --- a/src/lazyapplying.jl +++ b/src/lazyapplying.jl @@ -91,31 +91,3 @@ MemoryLayout(M::ApplyArray) = ApplyLayout(M.applied.f, MemoryLayout.(M.applied.a # _flatten(A::ApplyArray, B...) = _flatten(A.mul.args..., B...) # flatten(A::MulArray) = MulArray(Mul(_flatten(A.mul.args...))) - - - -const Add{Factors<:Tuple} = Applied{<:Any, typeof(+), Factors} - -size(M::Add, p::Int) = size(M)[p] -axes(M::Add, p::Int) = axes(M)[p] -ndims(M::Add) = ndims(first(M.args)) - -length(M::Add) = prod(size(M)) -size(M::Add) = length.(axes(M)) -axes(M::Add) = axes(first(M.args)) - - -eltype(M::Add) = Base._return_type(+, eltype.(M.args)) - -const AddArray{T,N,Factors<:Tuple} = ApplyArray{T,N,<:Add{Factors}} -const AddVector{T,Factors<:Tuple} = AddArray{T,1,Factors} -const AddMatrix{T,Factors<:Tuple} = AddArray{T,2,Factors} - -AddArray(factors...) = ApplyArray(+, factors...) - -""" - Add(A1, A2, …, AN) - -A lazy representation of `A1 + A2 + … + AN`; i.e., a shorthand for `applied(+, A1, A2, …, AN)`. -""" -Add(As...) = applied(+, As...) diff --git a/src/linalg/add.jl b/src/linalg/add.jl new file mode 100644 index 00000000..1a493a54 --- /dev/null +++ b/src/linalg/add.jl @@ -0,0 +1,48 @@ +#### +# These are special routines to make operations involving + +# more efficient +#### + + +const Add{Factors<:Tuple} = Applied{<:Any, typeof(+), Factors} + +size(M::Add, p::Int) = size(M)[p] +axes(M::Add, p::Int) = axes(M)[p] +ndims(M::Add) = ndims(first(M.args)) + +length(M::Add) = prod(size(M)) +size(M::Add) = length.(axes(M)) +axes(M::Add) = axes(first(M.args)) + + +eltype(M::Add) = Base._return_type(+, eltype.(M.args)) + +const AddArray{T,N,Factors<:Tuple} = ApplyArray{T,N,<:Add{Factors}} +const AddVector{T,Factors<:Tuple} = AddArray{T,1,Factors} +const AddMatrix{T,Factors<:Tuple} = AddArray{T,2,Factors} + +AddArray(factors...) = ApplyArray(+, factors...) + +""" + Add(A1, A2, …, AN) + +A lazy representation of `A1 + A2 + … + AN`; i.e., a shorthand for `applied(+, A1, A2, …, AN)`. +""" +Add(As...) = applied(+, As...) + + +for MulAdd_ in [MatMulMatAdd, MatMulVecAdd] + # `MulAdd{<:ApplyLayout{typeof(+)}}` cannot "win" against + # `MatMulMatAdd` and `MatMulVecAdd` hence `@eval`: + @eval function materialize!(M::$MulAdd_{<:ApplyLayout{typeof(+)}}) + α, A, B, β, C = M.α, M.A, M.B, M.β, M.C + if C ≡ B + B = copy(B) + end + lmul!(β, C) + for A in A.applied.args + C .= α .* Mul(A, B) .+ C + end + C + end +end diff --git a/src/linalg/blasmul.jl b/src/linalg/blasmul.jl index d2c57dad..27ae5916 100644 --- a/src/linalg/blasmul.jl +++ b/src/linalg/blasmul.jl @@ -202,22 +202,6 @@ function materialize!(M::MatMulVecAdd) default_blasmul!(α, A, B, iszero(β) ? false : β, C) end -for MulAdd_ in [MatMulMatAdd, MatMulVecAdd] - # `MulAdd{<:BroadcastLayout{typeof(+)}}` cannot "win" against - # `MatMulMatAdd` and `MatMulVecAdd` hence `@eval`: - @eval function materialize!(M::$MulAdd_{<:ApplyLayout{typeof(+)}}) - α, A, B, β, C = M.α, M.A, M.B, M.β, M.C - if C ≡ B - B = copy(B) - end - lmul!(β, C) - for A in A.applied.args - C .= α .* Mul(A, B) .+ C - end - C - end -end - # make copy to make sure always works @inline function _gemv!(tA, α, A, x, β, y) if x ≡ y diff --git a/src/linalg/linalg.jl b/src/linalg/linalg.jl index a9af6ef3..323c6724 100644 --- a/src/linalg/linalg.jl +++ b/src/linalg/linalg.jl @@ -3,3 +3,4 @@ include("lazymul.jl") include("blasbroadcasting.jl") include("blasmul.jl") include("inv.jl") +include("add.jl") From a47f0e5b6215bbcc2e11c62fad46a38d2d277351 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 24 Jan 2019 14:04:11 +0000 Subject: [PATCH 09/15] Support Mul(Add(...), b) --- src/lazyapplying.jl | 11 +++- src/linalg/add.jl | 5 ++ src/linalg/blasbroadcasting.jl | 1 - src/linalg/mul.jl | 34 +++++------- test/addtests.jl | 99 ++++++++++++++++++++++++++++++++++ test/multests.jl | 92 +------------------------------ test/runtests.jl | 1 + 7 files changed, 127 insertions(+), 116 deletions(-) create mode 100644 test/addtests.jl diff --git a/src/lazyapplying.jl b/src/lazyapplying.jl index 09fd8f01..8ff74791 100644 --- a/src/lazyapplying.jl +++ b/src/lazyapplying.jl @@ -1,7 +1,5 @@ - - abstract type ApplyStyle end struct DefaultApplyStyle <: ApplyStyle end struct LayoutApplyStyle{Layouts<:Tuple} <: ApplyStyle @@ -25,6 +23,15 @@ materialize(A::Applied{DefaultApplyStyle}) = A.f(materialize.(A.args)...) similar(M::Applied) = similar(M, eltype(M)) +struct ApplyBroadcastStyle <: BroadcastStyle end + +@inline copyto!(dest::AbstractArray, bc::Broadcasted{ApplyBroadcastStyle}) = + _copyto!(MemoryLayout(dest), dest, bc) +# Use default broacasting in general +@inline _copyto!(_, dest, bc::Broadcasted) = copyto!(dest, Broadcasted{Nothing}(bc.f, bc.args, bc.axes)) + +BroadcastStyle(::Type{<:Applied}) = ApplyBroadcastStyle() + struct MatrixFunctionStyle{F} <: ApplyStyle end diff --git a/src/linalg/add.jl b/src/linalg/add.jl index 1a493a54..94eac335 100644 --- a/src/linalg/add.jl +++ b/src/linalg/add.jl @@ -31,6 +31,11 @@ A lazy representation of `A1 + A2 + … + AN`; i.e., a shorthand for `applied(+, Add(As...) = applied(+, As...) +getindex(M::Add, k::Integer) = sum(getindex.(M.args, k)) +getindex(M::Add, k::Integer, j::Integer) = sum(getindex.(M.args, k, j)) +getindex(M::Add, k::CartesianIndex{1}) = M[convert(Int, k)] +getindex(M::Add, kj::CartesianIndex{2}) = M[kj[1], kj[2]] + for MulAdd_ in [MatMulMatAdd, MatMulVecAdd] # `MulAdd{<:ApplyLayout{typeof(+)}}` cannot "win" against # `MatMulMatAdd` and `MatMulVecAdd` hence `@eval`: diff --git a/src/linalg/blasbroadcasting.jl b/src/linalg/blasbroadcasting.jl index 954da474..b354f2d6 100644 --- a/src/linalg/blasbroadcasting.jl +++ b/src/linalg/blasbroadcasting.jl @@ -12,7 +12,6 @@ struct ArrayMulArrayStyle{StyleA, StyleB, p, q} <: BroadcastStyle end @inline copyto!(dest::AbstractArray, bc::Broadcasted{<:ArrayMulArrayStyle}) = _copyto!(MemoryLayout(dest), dest, bc) # Use default broacasting in general -@inline _copyto!(_, dest, bc::Broadcasted) = copyto!(dest, Broadcasted{Nothing}(bc.f, bc.args, bc.axes)) const BArrayMulArray{styleA, styleB, p, q, T, V} = Broadcasted{ArrayMulArrayStyle{styleA,styleB,p,q}, <:Any, typeof(identity), diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index 4d7bebcc..4e7a3430 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -108,16 +108,9 @@ gives an iterator containing the possible non-zero entries in the j-th column of colsupport(A, j) = colsupport(MemoryLayout(A), A, j) -function getindex(M::MatMulVec, k::Integer) - A,B = M.args - ret = zero(eltype(M)) - for j = rowsupport(A, k) - ret += A[k,j] * B[j] - end - ret -end -getindex(M::MatMulVec, k::CartesianIndex{1}) = M[convert(Int, k)] + + @@ -129,25 +122,18 @@ getindex(M::MatMulVec, k::CartesianIndex{1}) = M[convert(Int, k)] const MatMulMat{styleA, styleB, T, V} = ArrayMulArray{styleA, styleB, 2, 2, T, V} -function getindex(M::MatMulMat, k::Integer, j::Integer) - A,B = M.args - ret = zero(eltype(M)) - @inbounds for ℓ in (rowsupport(A,k) ∩ colsupport(B,j)) - ret += A[k,ℓ] * B[ℓ,j] - end - ret -end - -getindex(M::MatMulMat, kj::CartesianIndex{2}) = M[kj[1], kj[2]] #### # MulArray ##### +_mul(A) = A +_mul(A,B,C...) = Mul(A,B,C...) + function getindex(M::Mul, k::Integer) A,Bs = first(M.args), tail(M.args) - B = Mul(Bs) + B = _mul(Bs...) ret = zero(eltype(M)) for j = rowsupport(A, k) ret += A[k,j] * B[j] @@ -155,8 +141,12 @@ function getindex(M::Mul, k::Integer) ret end -_mul(A) = A -_mul(A,B,C...) = Mul(A,B,C...) + +getindex(M::Mul, k::CartesianIndex{1}) = M[convert(Int, k)] +getindex(M::Mul, kj::CartesianIndex{2}) = M[kj[1], kj[2]] + + + function getindex(M::Mul, k::Integer, j::Integer) A,Bs = first(M.args), tail(M.args) diff --git a/test/addtests.jl b/test/addtests.jl new file mode 100644 index 00000000..5c491a81 --- /dev/null +++ b/test/addtests.jl @@ -0,0 +1,99 @@ +using LazyArrays, Test +import LazyArrays: Add + + +@testset "Add" begin + @testset "Mul-Add" begin + A = Add(randn(5,5), randn(5,5)) + b = randn(5) + c = similar(b) + @test (c .= Mul(A, b)) ≈ A.args[1]*b + A.args[2]*b + end + + @testset "gemv Float64" begin + for A in (AddArray(randn(5,5), randn(5,5)), + AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))), + b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)) + + à = copy(A) + c = similar(b) + + c .= Mul(A,b) + @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c)) + + copyto!(c, Mul(A,b)) + @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c)) + + b̃ = copy(b) + copyto!(b̃, Mul(A,b̃)) + @test c ≈ b̃ + + c .= 2.0 .* Mul(A,b) + @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 0.0, similar(c)) + + c = copy(b) + c .= Mul(A,b) .+ c + @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 1.0, copy(b)) + + c = copy(b) + c .= Mul(A,b) .+ 2.0 .* c + @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 2.0, copy(b)) + + c = copy(b) + c .= 2.0 .* Mul(A,b) .+ c + @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 1.0, copy(b)) + + c = copy(b) + c .= 3.0 .* Mul(A,b) .+ 2.0 .* c + @test c ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, copy(b)) + + d = similar(c) + c = copy(b) + d .= 3.0 .* Mul(A,b) .+ 2.0 .* c + @test d ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, copy(b)) + end + end + + @testset "gemm" begin + for A in (AddArray(randn(5,5), randn(5,5)), + AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))), + B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), + view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)) + + à = copy(A) + C = similar(B) + + C .= Mul(A,B) + @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 0.0, similar(C)) + + B .= Mul(A,B) + @test C ≈ B + + C .= 2.0 .* Mul(A,B) + @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 0.0, similar(C)) + + C = copy(B) + C .= Mul(A,B) .+ C + @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 1.0, copy(B)) + + + C = copy(B) + C .= Mul(A,B) .+ 2.0 .* C + @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 2.0, copy(B)) + + C = copy(B) + C .= 2.0 .* Mul(A,B) .+ C + @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 1.0, copy(B)) + + + C = copy(B) + C .= 3.0 .* Mul(A,B) .+ 2.0 .* C + @test C ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, copy(B)) + + d = similar(C) + C = copy(B) + d .= 3.0 .* Mul(A,B) .+ 2.0 .* C + @test d ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, copy(B)) + end + end +end diff --git a/test/multests.jl b/test/multests.jl index 70429f75..46de0e86 100644 --- a/test/multests.jl +++ b/test/multests.jl @@ -686,96 +686,6 @@ import Base.Broadcast: materialize, materialize! M = Mul(A, Mul(B, C)) @test ndims(M) == ndims(typeof(M)) == 2 @test eltype(M) == Float64 - @test_broken all(copyto!(R1, M) .=== A*(B*C) .=== (R2 .= M)) - end -end - - -@testset "Add" begin - @testset "gemv Float64" begin - for A in (AddArray(randn(5,5), randn(5,5)), - AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))), - b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)) - - à = copy(A) - c = similar(b) - - c .= Mul(A,b) - @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c)) - - copyto!(c, Mul(A,b)) - @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c)) - - b̃ = copy(b) - copyto!(b̃, Mul(A,b̃)) - @test c ≈ b̃ - - c .= 2.0 .* Mul(A,b) - @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 0.0, similar(c)) - - c = copy(b) - c .= Mul(A,b) .+ c - @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 1.0, copy(b)) - - c = copy(b) - c .= Mul(A,b) .+ 2.0 .* c - @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 2.0, copy(b)) - - c = copy(b) - c .= 2.0 .* Mul(A,b) .+ c - @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 1.0, copy(b)) - - c = copy(b) - c .= 3.0 .* Mul(A,b) .+ 2.0 .* c - @test c ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, copy(b)) - - d = similar(c) - c = copy(b) - d .= 3.0 .* Mul(A,b) .+ 2.0 .* c - @test d ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, copy(b)) - end - end - - @testset "gemm" begin - for A in (AddArray(randn(5,5), randn(5,5)), - AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))), - B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), - view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)) - - à = copy(A) - C = similar(B) - - C .= Mul(A,B) - @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 0.0, similar(C)) - - B .= Mul(A,B) - @test C ≈ B - - C .= 2.0 .* Mul(A,B) - @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 0.0, similar(C)) - - C = copy(B) - C .= Mul(A,B) .+ C - @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 1.0, copy(B)) - - - C = copy(B) - C .= Mul(A,B) .+ 2.0 .* C - @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 2.0, copy(B)) - - C = copy(B) - C .= 2.0 .* Mul(A,B) .+ C - @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 1.0, copy(B)) - - - C = copy(B) - C .= 3.0 .* Mul(A,B) .+ 2.0 .* C - @test C ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, copy(B)) - - d = similar(C) - C = copy(B) - d .= 3.0 .* Mul(A,B) .+ 2.0 .* C - @test d ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, copy(B)) - end + @test all(copyto!(R1, M) .=== A*(B*C) .=== (R2 .= M)) end end diff --git a/test/runtests.jl b/test/runtests.jl index ac9ccdef..7f21e9a7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ include("memorylayouttests.jl") include("applytests.jl") include("multests.jl") include("ldivtests.jl") +include("addtests.jl") include("setoptests.jl") From 06c35b5103d63f813275fc7098d36d25d004b58a Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 24 Jan 2019 14:35:05 +0000 Subject: [PATCH 10/15] check dimensions for MulArray --- src/linalg/mul.jl | 6 +++++- test/multests.jl | 1 + 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index 4e7a3430..d7a5fdcc 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -164,6 +164,7 @@ const MulArray{T, N, MUL<:Mul} = ApplyArray{T, N, MUL} const MulVector{T, MUL<:Mul} = MulArray{T, 1, MUL} const MulMatrix{T, MUL<:Mul} = MulArray{T, 2, MUL} +const MulLayout{LAY} = ApplyLayout{typeof(*),LAY} MulArray{T,N}(M::MUL) where {T,N,MUL<:Mul} = MulArray{T,N,MUL}(M) MulArray{T}(M::Mul) where {T} = MulArray{T,ndims(M)}(M) @@ -171,7 +172,10 @@ MulArray(M::Mul) = MulArray{eltype(M)}(M) MulVector(M::Mul) = MulVector{eltype(M)}(M) MulMatrix(M::Mul) = MulMatrix{eltype(M)}(M) -MulArray(factors...) = MulArray(Mul(factors...)) +function MulArray(factors...) + checkdimensions(factors...) + MulArray(Mul(factors...)) +end MulArray{T}(factors...) where T = MulArray{T}(Mul(factors...)) MulArray{T,N}(factors...) where {T,N} = MulArray{T,N}(Mul(factors...)) MulVector(factors...) = MulVector(Mul(factors...)) diff --git a/test/multests.jl b/test/multests.jl index 46de0e86..07fd2580 100644 --- a/test/multests.jl +++ b/test/multests.jl @@ -669,6 +669,7 @@ import Base.Broadcast: materialize, materialize! A = randn(5,5) M = MulArray(A,A) @test Matrix(M) ≈ A^2 + @test_throws DimensionMismatch MulArray(randn(5,5), randn(4)) end @testset "#14" begin From 8a83714c5a7edc457aa4372a7be2e223bde77723 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 24 Jan 2019 16:11:30 +0000 Subject: [PATCH 11/15] Fix terterary Mul eltype --- src/LazyArrays.jl | 2 +- src/linalg/mul.jl | 10 +++++++++- test/multests.jl | 9 +++++++-- 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/src/LazyArrays.jl b/src/LazyArrays.jl index 911aaf1b..0926f659 100644 --- a/src/LazyArrays.jl +++ b/src/LazyArrays.jl @@ -28,7 +28,7 @@ import Base: ReinterpretArray, ReshapedArray, AbstractCartesianIndex, Slice, import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted, combine_eltypes, DefaultArrayStyle, instantiate, materialize, - materialize! + materialize!, eltypes import LinearAlgebra: AbstractTriangular, checksquare, pinv diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index d7a5fdcc..f61ed909 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -26,7 +26,15 @@ ndims(::Type{<:Mul{<:Any,Args}}) where Args = _mul_ndims(Args) length(M::Mul) = prod(size(M)) size(M::Mul) = length.(axes(M)) -eltype(M::Mul) = Base.promote_op(*, eltype.(M.args)...) + +@inline _mul_eltype(A) = A +@inline _mul_eltype(A, B) = Base.promote_op(*, A, B) +@inline _mul_eltype(A, B, C, D...) = _mul_eltype(Base.promote_op(*, A, B), C, D...) + +@inline _eltypes() = tuple() +@inline _eltypes(A, B...) = tuple(eltype(A), _eltypes(B...)...) + +@inline eltype(M::Mul) = _mul_eltype(_eltypes(M.args...)...) _mul_axes(ax1, ::Tuple{}) = (ax1,) _mul_axes(ax1, ::Tuple{<:Any}) = (ax1,) diff --git a/test/multests.jl b/test/multests.jl index 07fd2580..29312e87 100644 --- a/test/multests.jl +++ b/test/multests.jl @@ -3,10 +3,11 @@ import LazyArrays: MulAdd, MemoryLayout, DenseColumnMajor, DiagonalLayout, SymTr import Base.Broadcast: materialize, materialize! + @testset "Mul" begin @testset "eltype" begin - @inferred(eltype(Mul(zeros(Int,2,2), zeros(Float64,2)))) == Float64 - @inferred(eltype(Mul(zeros(ComplexF16,2,2),zeros(Int,2,2),zeros(Float64,2)))) == ComplexF64 + @test @inferred(eltype(Mul(zeros(Int,2,2), zeros(Float64,2)))) == Float64 + @test @inferred(eltype(Mul(zeros(ComplexF16,2,2),zeros(Int,2,2),zeros(Float64,2)))) == ComplexF64 v = Mul(zeros(Int,2,2), zeros(Float64,2)) A = Mul(zeros(Int,2,2), zeros(Float64,2,2)) @@ -18,6 +19,10 @@ import Base.Broadcast: materialize, materialize! à = materialize(A) @test à isa Matrix{Float64} @test à == zeros(2,2) + + A = randn(6,5); B = randn(5,5); C = randn(5,6) + M = Mul(A,B,C) + @test @inferred(eltype(M)) == Float64 end @testset "gemv Float64" begin From 51a9e333b0b8c255424cef2df36ec55f7a3bbc08 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 25 Jan 2019 09:08:19 +0000 Subject: [PATCH 12/15] fix inv --- src/linalg/inv.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linalg/inv.jl b/src/linalg/inv.jl index b6c06a2f..f03f997a 100644 --- a/src/linalg/inv.jl +++ b/src/linalg/inv.jl @@ -1,7 +1,7 @@ -const PInv{Style, Typ} = Applied{<:LayoutApplyStyle{Style}, typeof(pinv), Tuple{Typ}} -const Inv{Style, Typ} = Applied{<:LayoutApplyStyle{Style}, typeof(inv), Tuple{Typ}} +const PInv{Style, Typ} = Applied{<:LayoutApplyStyle{<:Tuple{Style}}, typeof(pinv), <:Tuple{Typ}} +const Inv{Style, Typ} = Applied{<:LayoutApplyStyle{<:Tuple{Style}}, typeof(inv), <:Tuple{Typ}} Inv(A) = applied(inv, A) PInv(A) = applied(pinv, A) From 9472a00478c961c6e2de6c68ed65719f3c0d6e04 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 28 Jan 2019 15:53:31 +0000 Subject: [PATCH 13/15] generalisations for ContinuumArrays --- src/lazyapplying.jl | 5 +++++ src/linalg/inv.jl | 7 ++++++- src/linalg/mul.jl | 9 ++++++--- 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/src/lazyapplying.jl b/src/lazyapplying.jl index 8ff74791..bf3ff362 100644 --- a/src/lazyapplying.jl +++ b/src/lazyapplying.jl @@ -6,6 +6,9 @@ struct LayoutApplyStyle{Layouts<:Tuple} <: ApplyStyle layouts::Layouts end +# Used for when a lazy version should be constructed on materialize +struct LazyArrayApplyStyle <: ApplyStyle end + ApplyStyle(f, args...) = DefaultApplyStyle() struct Applied{Style<:ApplyStyle, F, Args<:Tuple} @@ -79,6 +82,8 @@ IndexStyle(::ApplyArray{<:Any,1}) = IndexLinear() @propagate_inbounds getindex(A::ApplyArray{T,N}, kj::Vararg{Int,N}) where {T,N} = materialize(A.applied)[kj...] +materialize(A::Applied{LazyArrayApplyStyle}) = ApplyArray(A) + @inline copyto!(dest::AbstractArray, M::Applied) = _copyto!(MemoryLayout(dest), dest, M) @inline _copyto!(_, dest::AbstractArray, M::Applied) = copyto!(dest, materialize(M)) diff --git a/src/linalg/inv.jl b/src/linalg/inv.jl index f03f997a..94654883 100644 --- a/src/linalg/inv.jl +++ b/src/linalg/inv.jl @@ -23,6 +23,8 @@ end inv(A::Inv) = parent(A) pinv(A::Inv) = inv(A) +ndims(A::InvOrPInv) = ndims(parent(A)) + @@ -48,7 +50,10 @@ size(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) = (size(L.args[1], 2),) length(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) = size(L.args[1], 2) -eltype(M::Ldiv) = promote_type(eltype.(M.args)...) + +ndims(L::Applied{<:Any, typeof(\)}) = ndims(last(L.args)) +eltype(M::Applied{<:Any, typeof(\)}) = promote_type(Base.promote_op(inv, eltype(first(M.args))), + eltype(last(M.args))) struct ArrayLdivArrayStyle{StyleA, StyleB, p, q} <: BroadcastStyle end diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index f61ed909..98649c23 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -173,6 +173,7 @@ const MulVector{T, MUL<:Mul} = MulArray{T, 1, MUL} const MulMatrix{T, MUL<:Mul} = MulArray{T, 2, MUL} const MulLayout{LAY} = ApplyLayout{typeof(*),LAY} +MulLayout(layouts) = ApplyLayout(*, layouts) MulArray{T,N}(M::MUL) where {T,N,MUL<:Mul} = MulArray{T,N,MUL}(M) MulArray{T}(M::Mul) where {T} = MulArray{T,ndims(M)}(M) @@ -198,6 +199,8 @@ IndexStyle(::MulArray{<:Any,1}) = IndexLinear() adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.applied.args))...) transpose(A::MulArray) = MulArray(reverse(transpose.(A.applied.args))...) - -_flatten(A::MulArray, B...) = _flatten(A.mul.args..., B...) -flatten(A::MulArray) = MulArray(Mul(_flatten(A.mul.args...))) +_flatten() = () +_flatten(A, B...) = (A, _flatten(B...)...) +_flatten(A::Mul, B...) = _flatten(A.args..., B...) +_flatten(A::MulArray, B...) = _flatten(A.applied.args..., B...) +flatten(A::MulArray) = MulArray(Mul(_flatten(A.applied.args...)...)) From d39c0541db3d878aca6b0a427e2a670e24efd963 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 28 Jan 2019 16:02:39 +0000 Subject: [PATCH 14/15] Speed up compile time by stricter type definitions --- src/linalg/inv.jl | 6 +++--- src/linalg/mul.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/linalg/inv.jl b/src/linalg/inv.jl index 94654883..61ec487f 100644 --- a/src/linalg/inv.jl +++ b/src/linalg/inv.jl @@ -1,7 +1,7 @@ -const PInv{Style, Typ} = Applied{<:LayoutApplyStyle{<:Tuple{Style}}, typeof(pinv), <:Tuple{Typ}} -const Inv{Style, Typ} = Applied{<:LayoutApplyStyle{<:Tuple{Style}}, typeof(inv), <:Tuple{Typ}} +const PInv{Style, Typ} = Applied{LayoutApplyStyle{Tuple{Style}}, typeof(pinv), <:Tuple{Typ}} +const Inv{Style, Typ} = Applied{LayoutApplyStyle{Tuple{Style}}, typeof(inv), <:Tuple{Typ}} Inv(A) = applied(inv, A) PInv(A) = applied(pinv, A) @@ -37,7 +37,7 @@ eltype(A::InvOrPInv) = eltype(parent(A)) const Ldiv{StyleA, StyleB, AType, BType} = - Applied{<:LayoutApplyStyle{<:Tuple{StyleA, StyleB}}, typeof(\), <:Tuple{AType, BType}} + Applied{LayoutApplyStyle{Tuple{StyleA, StyleB}}, typeof(\), <:Tuple{AType, BType}} Ldiv(A, B) = applied(\, A, B) diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index 98649c23..d33a5355 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -7,7 +7,7 @@ function checkdimensions(A, B, C...) checkdimensions(B, C...) end -const Mul{Styles<:Tuple, Factors<:Tuple} = Applied{<:LayoutApplyStyle{Styles}, typeof(*), Factors} +const Mul{Styles<:Tuple, Factors<:Tuple} = Applied{LayoutApplyStyle{Styles}, typeof(*), Factors} ApplyStyle(::typeof(*), args::AbstractArray...) = LayoutApplyStyle(MemoryLayout.(args)) From 1c043508d0ad19d14f70c019d674b934846b0e28 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 10 Feb 2019 16:26:58 +0000 Subject: [PATCH 15/15] fix test, simplify materialize --- src/lazyapplying.jl | 14 +++++++------- src/linalg/mul.jl | 35 +---------------------------------- test/multests.jl | 3 ++- 3 files changed, 10 insertions(+), 42 deletions(-) diff --git a/src/lazyapplying.jl b/src/lazyapplying.jl index bf3ff362..bffe3135 100644 --- a/src/lazyapplying.jl +++ b/src/lazyapplying.jl @@ -18,9 +18,9 @@ struct Applied{Style<:ApplyStyle, F, Args<:Tuple} end applied(f, args...) = Applied(ApplyStyle(f, args...), f, args) -materialize(A::Applied{DefaultApplyStyle,<:Any,<:Tuple{<:Any}}) = - A.f(materialize(first(A.args))) -materialize(A::Applied{DefaultApplyStyle}) = A.f(materialize.(A.args)...) +_materialize(A::Applied, _) = A.f(materialize.(A.args)...) +materialize(M::Applied{<:LayoutApplyStyle}) = _materialize(M, axes(M)) +materialize(A::Applied) = A.f(materialize.(A.args)...) @@ -89,8 +89,8 @@ materialize(A::Applied{LazyArrayApplyStyle}) = ApplyArray(A) broadcastable(M::Applied) = M -# adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.mul.args))...) -# transpose(A::MulArray) = MulArray(reverse(transpose.(A.mul.args))...) +# adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.applied.args))...) +# transpose(A::MulArray) = MulArray(reverse(transpose.(A.applied.args))...) struct ApplyLayout{F, LAY} <: MemoryLayout @@ -101,5 +101,5 @@ end MemoryLayout(M::ApplyArray) = ApplyLayout(M.applied.f, MemoryLayout.(M.applied.args)) -# _flatten(A::ApplyArray, B...) = _flatten(A.mul.args..., B...) -# flatten(A::MulArray) = MulArray(Mul(_flatten(A.mul.args...))) +# _flatten(A::ApplyArray, B...) = _flatten(A.applied.args..., B...) +# flatten(A::MulArray) = MulArray(Mul(_flatten(A.applied.args...))) diff --git a/src/linalg/mul.jl b/src/linalg/mul.jl index d33a5355..70a7cc9c 100644 --- a/src/linalg/mul.jl +++ b/src/linalg/mul.jl @@ -43,30 +43,6 @@ axes(M::Mul) = _mul_axes(axes(first(M.args),1), axes(last(M.args))) axes(M::Mul{Tuple{}}) = () - - - -""" - lmaterialize(M::Mul) - -materializes arrays iteratively, left-to-right. -""" -lmaterialize(M::Mul) = _lmaterialize(M.args...) - -_lmaterialize(A, B) = materialize(Mul(A,B)) -_lmaterialize(A, B, C, D...) = _lmaterialize(materialize(Mul(A,B)), C, D...) - -""" - rmaterialize(M::Mul) - -materializes arrays iteratively, right-to-left. -""" -rmaterialize(M::Mul) = _rmaterialize(reverse(M.args)...) - -_rmaterialize(Z, Y) = materialize(Mul(Y,Z)) -_rmaterialize(Z, Y, X, W...) = _rmaterialize(materialize(Mul(Y,Z)), X, W...) - - # *(A::Mul, B::Mul) = materialize(Mul(A.args..., B.args...)) # *(A::Mul, B) = materialize(Mul(A.args..., B)) # *(A, B::Mul) = materialize(Mul(A, B.args...)) @@ -86,10 +62,7 @@ const ArrayMuls = Mul{<:Tuple, <:Tuple{Vararg{<:AbstractArray}}} similar(M::ArrayMuls, ::Type{T}, ::NTuple{N,OneTo{Int}}) where {T,N} = Array{T}(undef, size(M)) similar(M::ArrayMuls, ::Type{T}) where T = similar(M, T, axes(M)) _materialize(M::ArrayMulArray, _) = copyto!(similar(M), M) -_materialize(M::ArrayMuls, _) = lmaterialize(M) -_materialize(M::Mul, _) = lmaterialize(M) -_materialize(M::Mul2, _) = *(materialize.(M.args)...) -materialize(M::Mul) = _materialize(M, axes(M)) + #### @@ -198,9 +171,3 @@ IndexStyle(::MulArray{<:Any,1}) = IndexLinear() adjoint(A::MulArray) = MulArray(reverse(adjoint.(A.applied.args))...) transpose(A::MulArray) = MulArray(reverse(transpose.(A.applied.args))...) - -_flatten() = () -_flatten(A, B...) = (A, _flatten(B...)...) -_flatten(A::Mul, B...) = _flatten(A.args..., B...) -_flatten(A::MulArray, B...) = _flatten(A.applied.args..., B...) -flatten(A::MulArray) = MulArray(Mul(_flatten(A.applied.args...)...)) diff --git a/test/multests.jl b/test/multests.jl index 29312e87..3794f276 100644 --- a/test/multests.jl +++ b/test/multests.jl @@ -692,6 +692,7 @@ import Base.Broadcast: materialize, materialize! M = Mul(A, Mul(B, C)) @test ndims(M) == ndims(typeof(M)) == 2 @test eltype(M) == Float64 - @test all(copyto!(R1, M) .=== A*(B*C) .=== (R2 .= M)) + @test_skip all(copyto!(R1, M) .=== A*(B*C) .=== (R2 .= M)) + @test copyto!(R1, M) == A*(B*C) == (R2 .= M) end end