Move Out MemoryLayout, MulAdd, and L/Rmul to ArrayLayouts (#77)

* Only use tiled for finite arrays, better print for Vcat

* Move out code

* Move macros

* move _fill_lmul!

* Move zero!, scalarone/zero

* move factorizations

* Update LazyArrays.jl

* v0.14

Author: dlfivefifty

Project.toml

@@ -1,17 +1,20 @@
 name = "LazyArrays"
 uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02"
-version = "0.13.1"
+version = "0.14"
+
 
 [deps]
+ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-FillArrays = "0.7,0.8"
-MacroTools = "0.4.5,0.5"
-StaticArrays = "0.8,0.9,0.10,0.11,0.12"
+ArrayLayouts = "0.1"
+FillArrays = "0.8"
+MacroTools = "0.5"
+StaticArrays = "0.11, 0.12"
 julia = "1"
 
 [extras]

src/LazyArrays.jl

@@ -3,7 +3,7 @@ module LazyArrays
 # Use README as the docstring of the module:
 @doc read(joinpath(dirname(@__DIR__), "README.md"), String) LazyArrays
 
-using Base, Base.Broadcast, LinearAlgebra, FillArrays, StaticArrays
+using Base, Base.Broadcast, LinearAlgebra, FillArrays, StaticArrays, ArrayLayouts
 import LinearAlgebra.BLAS
 
 import Base: AbstractArray, AbstractMatrix, AbstractVector, 
@@ -45,6 +45,12 @@ import FillArrays: AbstractFill, getindex_value
 
 import StaticArrays: StaticArrayStyle
 
+import ArrayLayouts: MatMulVecAdd, MatMulMatAdd, MulAdd, Lmul, Rmul, Ldiv, 
+                        transposelayout, conjlayout, sublayout, triangularlayout, triangulardata,
+                        reshapedlayout, diagonallayout, adjointlayout,
+                        check_mul_axes, _mul_eltype, check_ldiv_axes, ldivaxes, colsupport, rowsupport,
+                        _fill_lmul!, @lazylmul, scalarone, scalarzero, fillzeros, zero!
+
 if VERSION < v"1.2-"
     import Base: has_offset_axes
     require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
@@ -56,7 +62,7 @@ export Mul, Applied, MulArray, MulVector, MulMatrix, InvMatrix, PInvMatrix,
         Hcat, Vcat, Kron, BroadcastArray, BroadcastMatrix, BroadcastVector, cache, Ldiv, Inv, PInv, Diff, Cumsum,
         applied, materialize, materialize!, ApplyArray, ApplyMatrix, ApplyVector, apply, ⋆, @~, LazyArray
 
-include("memorylayout.jl")
+
 include("lazyapplying.jl")
 include("lazybroadcasting.jl")
 include("linalg/linalg.jl")

src/lazyapplying.jl

@@ -68,7 +68,6 @@ similar(M::Applied) = similar(M, eltype(M))
 
 axes(A::Applied, j) = axes(A)[j]
 
-struct ApplyBroadcastStyle <: BroadcastStyle end
 struct ApplyArrayBroadcastStyle{N} <: Broadcast.AbstractArrayStyle{N} end
 ApplyArrayBroadcastStyle{N}(::Val{N}) where N = ApplyArrayBroadcastStyle{N}()
 
@@ -87,10 +86,7 @@ BroadcastStyle(::ApplyArrayBroadcastStyle{M}, b::DefaultArrayStyle{N}) where {M,
 similar(bc::Broadcasted{ApplyArrayBroadcastStyle{N}}, ::Type{ElType}) where {N,ElType} =
     similar(Array{ElType}, axes(bc))
 
-@inline function copyto!(dest::AbstractArray, bc::Broadcasted{ApplyBroadcastStyle}) 
-    @assert length(bc.args) == 1
-    copyto!(dest, first(bc.args))
-end
+
 @inline function copyto!(dest::AbstractArray, bc::Broadcasted{ApplyArrayBroadcastStyle{N}}) where N 
     if length(bc.args) == 1
         copyto!(dest, first(bc.args))
@@ -217,7 +213,7 @@ MemoryLayout(::Type{<:LazyArray}) = LazyArrayLayout()
 
 transposelayout(L::LazyLayout) = L
 conjlayout(L::LazyLayout) = L
-subarraylayout(L::LazyLayout, _) = L
+sublayout(L::LazyLayout, _) = L
 reshapedlayout(::LazyLayout, _) = LazyLayout()
 
 combine_mul_styles(::LazyLayout) = LazyArrayApplyStyle()
@@ -271,12 +267,6 @@ end
 # on the memory layout
 ###
 
-@inline sub_materialize(_, V) = Array(V)
-@inline sub_materialize(V::SubArray) = sub_materialize(MemoryLayout(typeof(V)), V)
-
-@inline lazy_getindex(A, I...) = sub_materialize(view(A, I...))
-
-
 @inline getindex(A::LazyMatrix, kr::Colon, jr::Colon) = lazy_getindex(A, kr, jr)
 @inline getindex(A::LazyMatrix, kr::Colon, jr::AbstractUnitRange) = lazy_getindex(A, kr, jr)
 @inline getindex(A::LazyMatrix, kr::AbstractUnitRange, jr::Colon) = lazy_getindex(A, kr, jr)

src/lazyconcat.jl

@@ -540,8 +540,8 @@ end
 # subarrays
 ###
 
-subarraylayout(::ApplyLayout{typeof(vcat)}, _) = ApplyLayout{typeof(vcat)}()
-subarraylayout(::ApplyLayout{typeof(hcat)}, _) = ApplyLayout{typeof(hcat)}()
+sublayout(::ApplyLayout{typeof(vcat)}, _) = ApplyLayout{typeof(vcat)}()
+sublayout(::ApplyLayout{typeof(hcat)}, _) = ApplyLayout{typeof(hcat)}()
 
 arguments(::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{<:Slice,<:Any}}) = 
     view.(arguments(parent(V)), Ref(:), Ref(parentindices(V)[2]))
@@ -623,13 +623,26 @@ materialize!(M::MatMulVecAdd{<:AbstractColumnMajor,<:ApplyLayout{typeof(vcat)}})
 
 ## print
 
+_replace_in_print_matrix(A::AbstractArray, k, j, s) = replace_in_print_matrix(A, k, j, s)
+_replace_in_print_matrix(_, k, j, s) = s
+
 function replace_in_print_matrix(f::Vcat{<:Any,1}, k::Integer, j::Integer, s::AbstractString)
     @assert j == 1
     κ = k
     for A in f.args
         n = length(A)
-        κ ≤ n && return replace_in_print_matrix(A, κ, 1, s)
+        κ ≤ n && return _replace_in_print_matrix(A, κ, 1, s)
         κ -= n
     end
     throw(BoundsError(f, k))
+end
+
+function replace_in_print_matrix(f::Vcat{<:Any,2}, k::Integer, j::Integer, s::AbstractString)
+    κ = k
+    for A in f.args
+        n = size(A,1)
+        κ ≤ n && return _replace_in_print_matrix(A, κ, j, s)
+        κ -= n
+    end
+    throw(BoundsError(f, (k,j)))
 end

src/linalg/add.jl

@@ -44,15 +44,6 @@ 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]]
 
-zero!(A::AbstractArray{T}) where T = fill!(A,zero(T))
-function zero!(A::AbstractArray{<:AbstractArray}) 
-    for a in A
-        zero!(a)
-    end
-    A
-end
-
-_fill_lmul!(β, A::AbstractArray{T}) where T = iszero(β) ? zero!(A) : lmul!(β, A)
 for MulAdd_ in [MatMulMatAdd, MatMulVecAdd]
     # `MulAdd{ApplyLayout{typeof(+)}}` cannot "win" against
     # `MatMulMatAdd` and `MatMulVecAdd` hence `@eval`:
@@ -89,7 +80,7 @@ end
 _view_tuple(a, b::Tuple) = view(a, b...)
 for op in (:+, :-)
     @eval begin
-        subarraylayout(a::ApplyLayout{typeof($op)}, _) = a
+        sublayout(a::ApplyLayout{typeof($op)}, _) = a
         arguments(::ApplyLayout{typeof($op)}, a::SubArray) =
             _view_tuple.(arguments(parent(a)), Ref(parentindices(a)))
         call(::ApplyLayout{typeof($op)}, a::SubArray) = $op            

src/linalg/diagonal.jl

@@ -1,5 +1,3 @@
-
-
 const PInv{Style, Typ} = Applied{Style, typeof(pinv), <:Tuple{Typ}}
 const Inv{Style, Typ} = Applied{Style, typeof(inv), <:Tuple{Typ}}
 
@@ -30,18 +28,6 @@ size(A::InvOrPInv, k) = size(A)[k]
 axes(A::InvOrPInv, k) = axes(A)[k]
 eltype(A::InvOrPInv) = Base.promote_op(inv, eltype(parent(A)))
 
-struct Ldiv{StyleA, StyleB, AType, BType}
-    A::AType
-    B::BType
-end
-
-Ldiv{StyleA, StyleB}(A::AType, B::BType) where {StyleA,StyleB,AType,BType} = 
-    Ldiv{StyleA,StyleB,AType,BType}(A,B)
-
-Ldiv(A::AType, B::BType) where {AType,BType} = 
-    Ldiv{typeof(MemoryLayout(AType)),typeof(MemoryLayout(BType)),AType,BType}(A, B)
-
-struct LdivBroadcastStyle <: BroadcastStyle end
 struct LdivApplyStyle <: ApplyStyle end
 
 
@@ -52,73 +38,16 @@ ldivapplystyle(_, ::LazyLayout) = LazyArrayApplyStyle()
 ApplyStyle(::typeof(\), ::Type{A}, ::Type{B}) where {A<:AbstractArray,B<:AbstractArray} = 
     ldivapplystyle(MemoryLayout(A), MemoryLayout(B))
 
-size(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractMatrix}) = (size(L.A, 2),size(L.B,2))
-size(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) = (size(L.A, 2),)
-axes(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractMatrix}) = (axes(L.A, 2),axes(L.B,2))
-axes(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) = (axes(L.A, 2),)    
-length(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) =size(L.A, 2)
-
-_ldivaxes(::Tuple{}, ::Tuple{}) = ()
-_ldivaxes(::Tuple{}, Bax::Tuple) = Bax
-_ldivaxes(::Tuple{<:Any}, ::Tuple{<:Any}) = ()
-_ldivaxes(::Tuple{<:Any}, Bax::Tuple{<:Any,<:Any}) = (OneTo(1),last(Bax))
-_ldivaxes(Aax::Tuple{<:Any,<:Any}, ::Tuple{<:Any}) = (last(Aax),)
-_ldivaxes(Aax::Tuple{<:Any,<:Any}, Bax::Tuple{<:Any,<:Any}) = (last(Aax),last(Bax))
-
-@inline ldivaxes(A, B) = _ldivaxes(axes(A), axes(B))
 
 axes(M::Applied{Style,typeof(\)}) where Style = ldivaxes(M.args...)
 axes(M::Applied{Style,typeof(\)}, p::Int)  where Style = axes(M)[p]
 size(M::Applied{Style,typeof(\)}) where Style = length.(axes(M))
 
-ndims(L::Ldiv) = ndims(last(L.args))
-eltype(M::Ldiv) = promote_type(Base.promote_op(inv, eltype(M.A)), eltype(M.B))
-
 @inline eltype(M::Applied{Style,typeof(\)}) where Style = eltype(Ldiv(M.args...))
 @inline ndims(M::Applied{Style,typeof(\)}) where Style = ndims(last(M.args))
 
-BroadcastStyle(::Type{<:Ldiv}) = ApplyBroadcastStyle()
-broadcastable(M::Ldiv) = M
-
-similar(A::Ldiv, ::Type{T}) where T = similar(Array{T}, axes(A))
-similar(A::Ldiv) = similar(A, eltype(A))
-
-function instantiate(L::Ldiv)
-    check_ldiv_axes(L.A, L.B)
-    Ldiv(instantiate(L.A), instantiate(L.B))
-end
-
-
-check_ldiv_axes(A, B) =
-    axes(A,1) == axes(B,1) || throw(DimensionMismatch("First axis of A, $(axes(A,1)), and first axis of B, $(axes(B,1)) must match"))
-
 check_applied_axes(A::Applied{<:Any,typeof(\)}) = check_ldiv_axes(A.args...)
 
-copy(M::Ldiv) = copyto!(similar(M), M)
-materialize(M::Ldiv) = copy(instantiate(M))
-
-
-
-_ldiv!(A, B) = ldiv!(factorize(A), B)
-_ldiv!(A::Factorization, B) = ldiv!(A, B)
-
-_ldiv!(dest, A, B) = ldiv!(dest, factorize(A), B)
-_ldiv!(dest, A::Factorization, B) = ldiv!(dest, A, B)
-
-
-materialize!(M::Ldiv) = _ldiv!(M.A, M.B)
-if VERSION ≥ v"1.1-pre"
-    copyto!(dest::AbstractArray, M::Ldiv) = _ldiv!(dest, M.A, M.B)
-else
-    copyto!(dest::AbstractArray, M::Ldiv) = _ldiv!(dest, M.A, copy(M.B))
-end
-
-const MatLdivVec{styleA, styleB, T, V} = Ldiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractVector{V}}
-const MatLdivMat{styleA, styleB, T, V} = Ldiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractMatrix{V}}
-const BlasMatLdivVec{styleA, styleB, T<:BlasFloat} = MatLdivVec{styleA, styleB, T, T}
-const BlasMatLdivMat{styleA, styleB, T<:BlasFloat} = MatLdivMat{styleA, styleB, T, T}
-
-
 ######
 # PInv/Inv
 ########

src/linalg/factorizations.jl

@@ -84,15 +84,6 @@ macro lazymul(Typ)
     esc(ret)
 end
 
-macro lazylmul(Typ)
-    esc(quote
-        LinearAlgebra.lmul!(A::$Typ, x::AbstractVector) = LazyArrays.materialize!(LazyArrays.Lmul(A,x))
-        LinearAlgebra.lmul!(A::$Typ, x::AbstractMatrix) = LazyArrays.materialize!(LazyArrays.Lmul(A,x))
-        LinearAlgebra.lmul!(A::$Typ, x::StridedVector) = LazyArrays.materialize!(LazyArrays.Lmul(A,x))
-        LinearAlgebra.lmul!(A::$Typ, x::StridedMatrix) = LazyArrays.materialize!(LazyArrays.Lmul(A,x))
-    end)
-end
-
 macro lazyldiv(Typ)
     esc(quote
         LinearAlgebra.ldiv!(A::$Typ, x::AbstractVector) = LazyArrays.materialize!(LazyArrays.Ldiv(A,x))

src/linalg/ldiv.jl renamed to src/linalg/inv.jl

@@ -1,10 +1,9 @@
 include("mul.jl")
-include("lazymul.jl")
 include("muladd.jl")
-include("lmul.jl")
-include("ldiv.jl")
+include("inv.jl")
+include("lazymul.jl")
 include("add.jl")
-include("factorizations.jl")
 
-include("diagonal.jl")
-include("triangular.jl")
+
+mulapplystyle(::TriangularLayout, ::AbstractStridedLayout) = LmulStyle()
+mulapplystyle(::AbstractStridedLayout, ::TriangularLayout) = RmulStyle()