Move over code from LazyArrays

Author: dlfivefifty

.travis.yml

@@ -0,0 +1,35 @@
+# Documentation: http://docs.travis-ci.com/user/languages/julia/
+language: julia
+os:
+  - linux
+  - osx
+  - windows
+julia:
+  - 1.0
+  - 1.1
+  - 1.2  
+  - 1.3  
+  - nightly
+matrix:
+  allow_failures:
+    - julia: nightly
+notifications:
+  email: false
+after_success:
+  - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())'
+
+# uncomment the following lines to override the default test script
+#script:
+#  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
+#  - julia -e 'Pkg.clone(pwd()); Pkg.build("MyPackage"); Pkg.test("MyPackage"; coverage=true)'
+
+# jobs:
+#   include:
+#     - stage: Documentation
+#       julia: 1.2
+#       script: julia --project=docs -e '
+#           using Pkg;
+#           Pkg.develop(PackageSpec(path=pwd()));
+#           Pkg.instantiate();
+#           include("docs/make.jl");'
+#       after_success: skip

Project.toml

@@ -0,0 +1,18 @@
+name = "ArrayLayouts"
+uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
+authors = ["Sheehan Olver <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[compat]
+FillArrays = "0.8"
+julia = "1"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]

src/ArrayLayouts.jl

@@ -0,0 +1,64 @@
+module ArrayLayouts
+using Base, Base.Broadcast, LinearAlgebra, FillArrays
+import LinearAlgebra.BLAS
+
+import Base: AbstractArray, AbstractMatrix, AbstractVector, 
+        ReinterpretArray, ReshapedArray, AbstractCartesianIndex, Slice,
+             RangeIndex, BroadcastStyle, copyto!, length, broadcastable, axes,
+             getindex, eltype, tail, IndexStyle, IndexLinear, getproperty,
+             *, +, -, /, \, ==, isinf, isfinite, sign, angle, show, isless,
+         fld, cld, div, min, max, minimum, maximum, mod,
+         <, ≤, >, ≥, promote_rule, convert, copy,
+         size, step, isempty, length, first, last, ndims,
+         getindex, setindex!, intersect, @_inline_meta, inv,
+         sort, sort!, issorted, sortperm, diff, cumsum, sum, in, broadcast,
+         eltype, parent, real, imag,
+         conj, transpose, adjoint, permutedims, vec,
+         exp, log, sqrt, cos, sin, tan, csc, sec, cot,
+                   cosh, sinh, tanh, csch, sech, coth,
+                   acos, asin, atan, acsc, asec, acot,
+                   acosh, asinh, atanh, acsch, asech, acoth, (:),
+         AbstractMatrix, AbstractArray, checkindex, unsafe_length, OneTo, one, zero,
+        to_shape, _sub2ind, print_matrix, print_matrix_row, print_matrix_vdots,
+      checkindex, Slice, @propagate_inbounds, @_propagate_inbounds_meta,
+      _in_range, _range, _rangestyle, Ordered,
+      ArithmeticWraps, floatrange, reverse, unitrange_last,
+      AbstractArray, AbstractVector, axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
+      similar, @_gc_preserve_end, @_gc_preserve_begin,
+      @nexprs, @ncall, @ntuple, tuple_type_tail,
+      all, any, isbitsunion, issubset, replace_in_print_matrix, replace_with_centered_mark
+
+import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted,
+                        combine_eltypes, DefaultArrayStyle, instantiate, materialize,
+                        materialize!, eltypes
+
+import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, Abuf, Bbuf, Cbuf, dot, factorize, qr, lu, cholesky, 
+                        norm2, norm1, normInf, normMinusInf
+
+import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex
+
+import FillArrays: AbstractFill, getindex_value
+
+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"))
+else
+    import Base: require_one_based_indexing    
+end     
+
+export materialize, materialize!, MulAdd, muladd!, Ldiv, Lmul, Rmul, MemoryLayout
+
+struct ApplyBroadcastStyle <: BroadcastStyle end
+@inline function copyto!(dest::AbstractArray, bc::Broadcasted{ApplyBroadcastStyle}) 
+    @assert length(bc.args) == 1
+    copyto!(dest, first(bc.args))
+end
+
+include("memorylayout.jl")
+include("muladd.jl")
+include("lmul.jl")
+include("ldiv.jl")
+include("diagonal.jl")
+include("triangular.jl")
+
+end

src/diagonal.jl

@@ -0,0 +1,31 @@
+
+rowsupport(::DiagonalLayout, _, k) = minimum(k):maximum(k)
+colsupport(::DiagonalLayout, _, j) = minimum(j):maximum(j)
+
+###
+# Lmul
+####
+
+# Diagonal multiplication never changes structure
+similar(M::Lmul{<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.B, T, axes)
+# equivalent to rescaling
+function materialize!(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}})
+    M.B .= getindex_value(M.A.diag) .* M.B
+    M.B
+end
+
+copy(M::Lmul{<:DiagonalLayout,<:DiagonalLayout}) = Diagonal(M.A.diag .* M.B.diag)
+copy(M::Lmul{<:DiagonalLayout}) = M.A.diag .* M.B
+copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}}) = getindex_value(M.A.diag) .* M.B
+copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = Diagonal(getindex_value(M.A.diag) .* M.B.diag)
+
+# Diagonal multiplication never changes structure
+similar(M::Rmul{<:Any,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.A, T, axes)
+# equivalent to rescaling
+function materialize!(M::Rmul{<:Any,<:DiagonalLayout{<:AbstractFillLayout}})
+    M.A .= M.A .* getindex_value(M.B.diag)
+    M.A
+end
+
+copy(M::Rmul{<:Any,<:DiagonalLayout}) = M.A .* permutedims(M.B.diag)
+copy(M::Rmul{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) =  M.A .* getindex_value(M.B.diag)

src/ldiv.jl

@@ -0,0 +1,73 @@
+
+
+
+
+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
+
+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))
+
+ndims(L::Ldiv) = ndims(last(L.args))
+eltype(M::Ldiv) = promote_type(Base.promote_op(inv, eltype(M.A)), eltype(M.B))
+
+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"))
+
+
+
+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}

src/lmul.jl

@@ -0,0 +1,60 @@
+for Typ in (:Lmul, :Rmul)
+    @eval begin
+        struct $Typ{StyleA, StyleB, TypeA, TypeB}
+            A::TypeA
+            B::TypeB
+        end
+
+        $Typ(A::TypeA, B::TypeB) where {TypeA,TypeB} = $Typ{typeof(MemoryLayout(TypeA)),typeof(MemoryLayout(TypeB)),TypeA,TypeB}(A,B)
+
+        BroadcastStyle(::Type{<:$Typ}) = ApplyBroadcastStyle()
+        broadcastable(M::$Typ) = M
+
+        eltype(::$Typ{<:Any,<:Any,A,B}) where {A,B} = promote_type(eltype(A), eltype(B))
+        size(M::$Typ, p::Int) = size(M)[p]
+        axes(M::$Typ, p::Int) = axes(M)[p]
+        length(M::$Typ) = prod(size(M))
+        size(M::$Typ) = map(length,axes(M))
+        axes(M::$Typ) = (axes(M.A,1),axes(M.B,2))
+
+        similar(M::$Typ, ::Type{T}, axes) where {T,N} = similar(Array{T}, axes)
+        similar(M::$Typ, ::Type{T}) where T = similar(M, T, axes(M))
+        similar(M::$Typ) = similar(M, eltype(M))
+    end
+end
+
+
+
+const MatLmulVec{StyleA,StyleB} = Lmul{StyleA,StyleB,<:AbstractMatrix,<:AbstractVector}
+const MatLmulMat{StyleA,StyleB} = Lmul{StyleA,StyleB,<:AbstractMatrix,<:AbstractMatrix}
+
+const BlasMatLmulVec{StyleA,StyleB,T<:BlasFloat} = Lmul{StyleA,StyleB,<:AbstractMatrix{T},<:AbstractVector{T}}
+const BlasMatLmulMat{StyleA,StyleB,T<:BlasFloat} = Lmul{StyleA,StyleB,<:AbstractMatrix{T},<:AbstractMatrix{T}}
+
+const MatRmulMat{StyleA,StyleB} = Rmul{StyleA,StyleB,<:AbstractMatrix,<:AbstractMatrix}
+const BlasMatRmulMat{StyleA,StyleB,T<:BlasFloat} = Rmul{StyleA,StyleB,<:AbstractMatrix{T},<:AbstractMatrix{T}}
+
+
+####
+# LMul materialize
+####
+
+axes(M::MatLmulVec) = (axes(M.A,1),)
+
+
+copy(M::Lmul) = materialize!(Lmul(M.A,copy(M.B)))
+copy(M::Rmul) = materialize!(Rmul(copy(M.A),M.B))
+
+@inline function copyto!(dest::AbstractArray, M::Lmul)
+    M.B ≡ dest || copyto!(dest, M.B)
+    materialize!(Lmul(M.A,dest))
+end
+
+@inline function copyto!(dest::AbstractArray, M::Rmul)
+    M.A ≡ dest || copyto!(dest, M.A)
+    materialize!(Rmul(dest,M.B))
+end
+
+materialize!(M::Lmul) = lmul!(M.A,M.B)
+materialize!(M::Rmul) = rmul!(M.A,M.B)
+