From e1f66e269a671c7b6c5e6fcd77f8619681a88c46 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 4 Aug 2023 18:11:11 +0200 Subject: [PATCH 1/7] WIP: make FillArrays.jl a weakdep --- Project.toml | 11 ++- ext/ArrayLayoutsFillArraysExt.jl | 113 +++++++++++++++++++++++++++++++ src/ArrayLayouts.jl | 27 +++++--- src/diagonal.jl | 26 ------- src/memorylayout.jl | 11 --- src/mul.jl | 11 --- src/muladd.jl | 18 ----- 7 files changed, 138 insertions(+), 79 deletions(-) create mode 100644 ext/ArrayLayoutsFillArraysExt.jl diff --git a/Project.toml b/Project.toml index d1fc06a..ee851ba 100644 --- a/Project.toml +++ b/Project.toml @@ -1,22 +1,29 @@ name = "ArrayLayouts" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" authors = ["Sheehan Olver "] -version = "1.1.1" +version = "1.2.0" [deps] FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +[weakdeps] +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" + +[extensions] +ArrayLayoutsFillArraysExt = "FillArrays" + [compat] FillArrays = "1.2.1" julia = "1.6" [extras] Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Base64", "Random", "StableRNGs", "Test"] +test = ["Base64", "FillArrays", "Random", "StableRNGs", "Test"] diff --git a/ext/ArrayLayoutsFillArraysExt.jl b/ext/ArrayLayoutsFillArraysExt.jl new file mode 100644 index 0000000..14e16be --- /dev/null +++ b/ext/ArrayLayoutsFillArraysExt.jl @@ -0,0 +1,113 @@ +module ArrayLayoutsFillArraysExt + +using FillArrays +using FillArrays: AbstractFill, getindex_value + +using ArrayLayouts +using ArrayLayouts: OnesLayout, Mul, MulAdd +import ArrayLayouts: MemoryLayout, _copyto!, sub_materialize, diagonaldata +export layoutfillmul, mulzeros + +import Base: copy, * +import Base.Broadcast: materialize! +import LinearAlgebra +using LinearAlgebra: Adjoint, Transpose, Symmetric, Hermitian, Diagonal, + AdjointAbsVec, TransposeAbsVec, UniformScaling + +macro layoutfillmul(Typ) + ret = quote + (*)(A::LinearAlgebra.AdjointAbsVec{<:Any,<:Zeros{<:Any,1}}, B::$Typ) = ArrayLayouts.mul(A,B) + (*)(A::LinearAlgebra.TransposeAbsVec{<:Any,<:Zeros{<:Any,1}}, B::$Typ) = ArrayLayouts.mul(A,B) + (*)(A::LinearAlgebra.Transpose{T,<:$Typ}, B::Zeros{T,1}) where T<:Real = ArrayLayouts.mul(A,B) + end + for Mod in (:Adjoint, :Transpose, :Symmetric, :Hermitian) + ret = quote + $ret + + (*)(A::$Mod{<:Any,<:$Typ}, B::Zeros{<:Any,1}) = ArrayLayouts.mul(A,B) + end + end + esc(ret) +end + +@layoutfillmul LayoutMatrix + +*(a::Zeros{<:Any,2}, b::LayoutMatrix) = FillArrays.mult_zeros(a, b) +*(a::LayoutMatrix, b::Zeros{<:Any,2}) = FillArrays.mult_zeros(a, b) +*(a::LayoutMatrix, b::Zeros{<:Any,1}) = FillArrays.mult_zeros(a, b) +*(a::Transpose{T, <:LayoutMatrix{T}} where T, b::Zeros{<:Any, 2}) = FillArrays.mult_zeros(a, b) +*(a::Adjoint{T, <:LayoutMatrix{T}} where T, b::Zeros{<:Any, 2}) = FillArrays.mult_zeros(a, b) +*(A::Adjoint{<:Any, <:Zeros{<:Any,1}}, B::Diagonal{<:Any,<:LayoutVector}) = (B' * A')' +*(A::Transpose{<:Any, <:Zeros{<:Any,1}}, B::Diagonal{<:Any,<:LayoutVector}) = transpose(transpose(B) * transpose(A)) + +# equivalent to rescaling +function materialize!(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}}) + M.B .= getindex_value(M.A.diag) .* M.B + M.B +end +# equivalent to rescaling +function materialize!(M::Rmul{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) + M.A .= M.A .* getindex_value(M.B.diag) + M.A +end + +copy(M::Ldiv{<:DiagonalLayout{<:AbstractFillLayout}}) = inv(getindex_value(M.A.diag)) .* M.B +copy(M::Ldiv{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = diagonal(inv(getindex_value(M.A.diag)) .* M.B.diag) + +copy(M::Rdiv{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A .* inv(getindex_value(M.B.diag)) +copy(M::Rdiv{<:DiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = diagonal(M.A.diag .* inv(getindex_value(M.B.diag))) + +copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}}) = getindex_value(diagonaldata(M.A)) * M.B +copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B +copy(M::Rmul{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) +copy(M::Rmul{<:DualLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) + +copy(M::Rmul{<:BidiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) +copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:BidiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B +copy(M::Rmul{<:TridiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) +copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:TridiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B +copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) +copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:SymTridiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B + +MemoryLayout(::Type{<:AbstractFill}) = FillLayout() +MemoryLayout(::Type{<:Zeros}) = ZerosLayout() +MemoryLayout(::Type{<:Ones}) = OnesLayout() + +_copyto!(_, ::AbstractFillLayout, dest::AbstractArray{<:Any,N}, src::AbstractArray{<:Any,N}) where N = + fill!(dest, getindex_value(src)) + + _fill_copyto!(dest, C::Zeros) = zero!(dest) # exploit special fill! overload + +sub_materialize(::AbstractFillLayout, V, ax) = Fill(getindex_value(V), ax) +sub_materialize(::ZerosLayout, V, ax) = Zeros{eltype(V)}(ax) +sub_materialize(::OnesLayout, V, ax) = Ones{eltype(V)}(ax) + +*(x::AdjointAbsVec{<:Any,<:Zeros{<:Any,1}}, D::Diagonal, y::LayoutVector) = FillArrays._triple_zeromul(x, D, y) +*(x::TransposeAbsVec{<:Any,<:Zeros{<:Any,1}}, D::Diagonal, y::LayoutVector) = FillArrays._triple_zeromul(x, D, y) + +@inline LinearAlgebra.dot(a::LayoutVector, b::AbstractFill{<:Any,1}) = FillArrays._fill_dot_rev(a,b) +@inline LinearAlgebra.dot(a::AbstractFill{<:Any,1}, b::LayoutVector) = FillArrays._fill_dot(a,b) + +# TODO: Remove unnecessary _apply +_apply(_, _, op, Λ::UniformScaling, A::AbstractMatrix) = op(Diagonal(Fill(Λ.λ,(axes(A,1),))), A) +_apply(_, _, op, A::AbstractMatrix, Λ::UniformScaling) = op(A, Diagonal(Fill(Λ.λ,(axes(A,1),)))) + +# equivalent to rescaling +function materialize!(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout}}) + checkdimensions(M) + M.C .= (M.α * getindex_value(M.A.diag)) .* M.B .+ M.β .* M.C + M.C +end + +function materialize!(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) + checkdimensions(M) + M.C .= M.α .* M.A .* getindex_value(M.B.diag) .+ M.β .* M.C + M.C +end + +fillzeros(::Type{T}, ax) where T<:Number = Zeros{T}(ax) +mulzeros(::Type{T}, M) where T<:Number = fillzeros(T, axes(M)) +mulzeros(::Type{T}, M::Mul{<:DualLayout,<:Any,<:Adjoint}) where T<:Number = fillzeros(T, axes(M,2))' +mulzeros(::Type{T}, M::Mul{<:DualLayout,<:Any,<:Transpose}) where T<:Number = transpose(fillzeros(T, axes(M,2))) + +end diff --git a/src/ArrayLayouts.jl b/src/ArrayLayouts.jl index cfc820c..0aa4a08 100644 --- a/src/ArrayLayouts.jl +++ b/src/ArrayLayouts.jl @@ -1,6 +1,6 @@ module ArrayLayouts using Base: _typed_hcat -using Base, Base.Broadcast, LinearAlgebra, FillArrays, SparseArrays +using Base, Base.Broadcast, LinearAlgebra, SparseArrays using LinearAlgebra.BLAS using Base: AbstractCartesianIndex, OneTo, oneto, RangeIndex, ReinterpretArray, ReshapedArray, @@ -25,8 +25,6 @@ using LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex AdjointQtype{T} = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ{T} : Adjoint{T,<:AbstractQ} -using FillArrays: AbstractFill, getindex_value, axes_print_matrix_row - using Base: require_one_based_indexing export materialize, materialize!, MulAdd, muladd!, Ldiv, Rdiv, Lmul, Rmul, Dot, @@ -121,6 +119,11 @@ include("diagonal.jl") include("triangular.jl") include("factorizations.jl") +@static if !isdefined(Base, :get_extension) + include("../ext/ArrayLayoutsFillArraysExt.jl") +end + + # Extend this function if you're only looking to dispatch on the axes @inline sub_materialize_axes(V, _) = Array(V) @inline sub_materialize(_, V, ax) = sub_materialize_axes(V, ax) @@ -197,14 +200,6 @@ getindex(A::LayoutVector, kr::Colon) = layout_getindex(A, kr) getindex(A::AdjOrTrans{<:Any,<:LayoutVector}, kr::Integer, jr::Colon) = layout_getindex(A, kr, jr) getindex(A::AdjOrTrans{<:Any,<:LayoutVector}, kr::Integer, jr::AbstractVector) = layout_getindex(A, kr, jr) -*(a::Zeros{<:Any,2}, b::LayoutMatrix) = FillArrays.mult_zeros(a, b) -*(a::LayoutMatrix, b::Zeros{<:Any,2}) = FillArrays.mult_zeros(a, b) -*(a::LayoutMatrix, b::Zeros{<:Any,1}) = FillArrays.mult_zeros(a, b) -*(a::Transpose{T, <:LayoutMatrix{T}} where T, b::Zeros{<:Any, 2}) = FillArrays.mult_zeros(a, b) -*(a::Adjoint{T, <:LayoutMatrix{T}} where T, b::Zeros{<:Any, 2}) = FillArrays.mult_zeros(a, b) -*(A::Adjoint{<:Any, <:Zeros{<:Any,1}}, B::Diagonal{<:Any,<:LayoutVector}) = (B' * A')' -*(A::Transpose{<:Any, <:Zeros{<:Any,1}}, B::Diagonal{<:Any,<:LayoutVector}) = transpose(transpose(B) * transpose(A)) - *(A::Diagonal{<:Any,<:LayoutVector}, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B) *(A::Diagonal{<:Any,<:LayoutVector}, B::AbstractMatrix) = mul(A, B) *(A::AbstractMatrix, B::Diagonal{<:Any,<:LayoutVector}) = mul(A, B) @@ -378,6 +373,16 @@ Base.replace_in_print_matrix(A::Union{LayoutVector, UnitLowerTriangular{<:Any,<:AdjOrTrans{<:Any,<:LayoutMatrix}}}, i::Integer, j::Integer, s::AbstractString) = layout_replace_in_print_matrix(MemoryLayout(A), A, i, j, s) +if VERSION < v"1.8-" + axes_print_matrix_row(lay, io, X, A, i, cols, sep) = + Base.invoke(Base.print_matrix_row, Tuple{IO,AbstractVecOrMat,Vector,Integer,AbstractVector,AbstractString}, + io, X, A, i, cols, sep) +else + axes_print_matrix_row(lay, io, X, A, i, cols, sep, idxlast::Integer=last(axes(X, 2))) = + Base.invoke(Base.print_matrix_row, Tuple{IO,AbstractVecOrMat,Vector,Integer,AbstractVector,AbstractString,Integer}, + io, X, A, i, cols, sep, idxlast) +end + Base.print_matrix_row(io::IO, X::Union{LayoutMatrix, LayoutVector, diff --git a/src/diagonal.jl b/src/diagonal.jl index 1a3e26d..22b86b0 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -8,11 +8,6 @@ mulreduce(M::Mul{<:Any,<:DiagonalLayout}) = Rmul(M) # 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(diagonaldata(M.A) .* diagonaldata(M.B)) @@ -32,11 +27,6 @@ end # 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 function materialize!(M::Ldiv{<:DiagonalLayout}) @@ -46,13 +36,9 @@ end copy(M::Ldiv{<:DiagonalLayout,<:DiagonalLayout}) = diagonal(M.A.diag .\ M.B.diag) copy(M::Ldiv{<:DiagonalLayout}) = M.A.diag .\ M.B -copy(M::Ldiv{<:DiagonalLayout{<:AbstractFillLayout}}) = inv(getindex_value(M.A.diag)) .* M.B -copy(M::Ldiv{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = diagonal(inv(getindex_value(M.A.diag)) .* M.B.diag) copy(M::Rdiv{<:DiagonalLayout,<:DiagonalLayout}) = diagonal(M.A.diag .* inv.(M.B.diag)) copy(M::Rdiv{<:Any,<:DiagonalLayout}) = M.A .* inv.(permutedims(M.B.diag)) -copy(M::Rdiv{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A .* inv(getindex_value(M.B.diag)) -copy(M::Rdiv{<:DiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = diagonal(M.A.diag .* inv(getindex_value(M.B.diag))) @@ -71,18 +57,6 @@ copy(M::Lmul{DiagonalLayout{OnesLayout},DiagonalLayout{OnesLayout}}) = _copy_oft copy(M::Rmul{<:Any,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M)) copy(M::Rmul{<:DualLayout,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M)) -copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}}) = getindex_value(diagonaldata(M.A)) * M.B -copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B -copy(M::Rmul{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) -copy(M::Rmul{<:DualLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) - -copy(M::Rmul{<:BidiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) -copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:BidiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B -copy(M::Rmul{<:TridiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) -copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:TridiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B -copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B)) -copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:SymTridiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B - copy(M::Rmul{<:BidiagonalLayout,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M)) copy(M::Lmul{DiagonalLayout{OnesLayout},<:BidiagonalLayout}) = _copy_oftype(M.B, eltype(M)) diff --git a/src/memorylayout.jl b/src/memorylayout.jl index 989e5f1..7f797f7 100644 --- a/src/memorylayout.jl +++ b/src/memorylayout.jl @@ -531,23 +531,12 @@ struct ZerosLayout <: AbstractFillLayout end struct OnesLayout <: AbstractFillLayout end struct EyeLayout <: MemoryLayout end -MemoryLayout(::Type{<:AbstractFill}) = FillLayout() -MemoryLayout(::Type{<:Zeros}) = ZerosLayout() -MemoryLayout(::Type{<:Ones}) = OnesLayout() - # all sub arrays are same sublayout(L::AbstractFillLayout, inds::Type) = L reshapedlayout(L::AbstractFillLayout, _) = L adjointlayout(::Type, L::AbstractFillLayout) = L transposelayout(L::AbstractFillLayout) = L -_copyto!(_, ::AbstractFillLayout, dest::AbstractArray{<:Any,N}, src::AbstractArray{<:Any,N}) where N = - fill!(dest, getindex_value(src)) - -sub_materialize(::AbstractFillLayout, V, ax) = Fill(getindex_value(V), ax) -sub_materialize(::ZerosLayout, V, ax) = Zeros{eltype(V)}(ax) -sub_materialize(::OnesLayout, V, ax) = Ones{eltype(V)}(ax) - abstract type AbstractBandedLayout <: MemoryLayout end abstract type AbstractTridiagonalLayout <: AbstractBandedLayout end diff --git a/src/mul.jl b/src/mul.jl index 6548464..6c1f749 100644 --- a/src/mul.jl +++ b/src/mul.jl @@ -213,9 +213,6 @@ macro layoutmul(Typ) (*)(A::AbstractMatrix, B::$Typ) = ArrayLayouts.mul(A,B) (*)(A::LinearAlgebra.AdjointAbsVec, B::$Typ) = ArrayLayouts.mul(A,B) (*)(A::LinearAlgebra.TransposeAbsVec, B::$Typ) = ArrayLayouts.mul(A,B) - (*)(A::LinearAlgebra.AdjointAbsVec{<:Any,<:Zeros{<:Any,1}}, B::$Typ) = ArrayLayouts.mul(A,B) - (*)(A::LinearAlgebra.TransposeAbsVec{<:Any,<:Zeros{<:Any,1}}, B::$Typ) = ArrayLayouts.mul(A,B) - (*)(A::LinearAlgebra.Transpose{T,<:$Typ}, B::Zeros{T,1}) where T<:Real = ArrayLayouts.mul(A,B) (*)(A::LinearAlgebra.AbstractQ, B::$Typ) = ArrayLayouts.mul(A,B) (*)(A::$Typ, B::LinearAlgebra.AbstractQ) = ArrayLayouts.mul(A,B) @@ -278,7 +275,6 @@ macro layoutmul(Typ) (*)(A::LinearAlgebra.TransposeAbsVec, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B) (*)(A::$Mod{<:Any,<:$Typ}, B::AbstractVector) = ArrayLayouts.mul(A,B) (*)(A::$Mod{<:Any,<:$Typ}, B::ArrayLayouts.LayoutVector) = ArrayLayouts.mul(A,B) - (*)(A::$Mod{<:Any,<:$Typ}, B::Zeros{<:Any,1}) = ArrayLayouts.mul(A,B) (*)(A::$Mod{<:Any,<:$Typ}, B::$Typ) = ArrayLayouts.mul(A,B) (*)(A::$Typ, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B) @@ -306,8 +302,6 @@ end *(x::Transpose{<:Any,<:LayoutVector}, D::Diagonal{<:Any,<:LayoutVector}) = mul(x, D) *(x::AdjointAbsVec, D::Diagonal, y::LayoutVector) = x * mul(D,y) *(x::TransposeAbsVec, D::Diagonal, y::LayoutVector) = x * mul(D,y) -*(x::AdjointAbsVec{<:Any,<:Zeros{<:Any,1}}, D::Diagonal, y::LayoutVector) = FillArrays._triple_zeromul(x, D, y) -*(x::TransposeAbsVec{<:Any,<:Zeros{<:Any,1}}, D::Diagonal, y::LayoutVector) = FillArrays._triple_zeromul(x, D, y) *(A::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}, B::UpperOrLowerTriangular{<:Any,<:LayoutMatrix}) = mul(A, B) @@ -348,8 +342,6 @@ dot(a, b) = materialize(Dot(a, b)) @inline LinearAlgebra.dot(a::LayoutArray, b::LayoutArray) = dot(a,b) @inline LinearAlgebra.dot(a::LayoutArray, b::AbstractArray) = dot(a,b) @inline LinearAlgebra.dot(a::AbstractArray, b::LayoutArray) = dot(a,b) -@inline LinearAlgebra.dot(a::LayoutVector, b::AbstractFill{<:Any,1}) = FillArrays._fill_dot_rev(a,b) -@inline LinearAlgebra.dot(a::AbstractFill{<:Any,1}, b::LayoutVector) = FillArrays._fill_dot(a,b) @inline LinearAlgebra.dot(a::LayoutArray{<:Number}, b::SparseArrays.SparseVectorUnion{<:Number}) = dot(a,b) @inline LinearAlgebra.dot(a::SparseArrays.SparseVectorUnion{<:Number}, b::LayoutArray{<:Number}) = dot(a,b) @@ -369,9 +361,6 @@ LinearAlgebra.dot(x::AbstractVector, A::Symmetric{<:Real,<:LayoutMatrix}, y::Abs # allow overloading for infinite or lazy case @inline _power_by_squaring(_, _, A, p) = invoke(Base.power_by_squaring, Tuple{AbstractMatrix,Integer}, A, p) -# TODO: Remove unnecessary _apply -_apply(_, _, op, Λ::UniformScaling, A::AbstractMatrix) = op(Diagonal(Fill(Λ.λ,(axes(A,1),))), A) -_apply(_, _, op, A::AbstractMatrix, Λ::UniformScaling) = op(A, Diagonal(Fill(Λ.λ,(axes(A,1),)))) for Typ in (:LayoutMatrix, :(Symmetric{<:Any,<:LayoutMatrix}), :(Hermitian{<:Any,<:LayoutMatrix}), :(Adjoint{<:Any,<:LayoutMatrix}), :(Transpose{<:Any,<:LayoutMatrix})) diff --git a/src/muladd.jl b/src/muladd.jl index 6087643..dff671f 100644 --- a/src/muladd.jl +++ b/src/muladd.jl @@ -72,7 +72,6 @@ materialize(M::MulAdd) = copy(instantiate(M)) copy(M::MulAdd) = copyto!(similar(M), M) _fill_copyto!(dest, C) = copyto!(dest, C) -_fill_copyto!(dest, C::Zeros) = zero!(dest) # exploit special fill! overload @inline copyto!(dest::AbstractArray{T}, M::MulAdd) where T = muladd!(M.α, unalias(dest,M.A), unalias(dest,M.B), M.β, _fill_copyto!(dest, M.C)) @@ -362,18 +361,6 @@ materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:AbstractRowMajor},<:Abstrac similar(M::MulAdd{<:DiagonalLayout,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.B, T, axes) similar(M::MulAdd{<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.B, T, axes) similar(M::MulAdd{<:Any,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.A, T, axes) -# equivalent to rescaling -function materialize!(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout}}) - checkdimensions(M) - M.C .= (M.α * getindex_value(M.A.diag)) .* M.B .+ M.β .* M.C - M.C -end - -function materialize!(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) - checkdimensions(M) - M.C .= M.α .* M.A .* getindex_value(M.B.diag) .+ M.β .* M.C - M.C -end BroadcastStyle(::Type{<:MulAdd}) = ApplyBroadcastStyle() @@ -383,11 +370,6 @@ scalarone(::Type{A}) where {A<:AbstractArray} = scalarone(eltype(A)) scalarzero(::Type{T}) where T = zero(T) scalarzero(::Type{A}) where {A<:AbstractArray} = scalarzero(eltype(A)) -fillzeros(::Type{T}, ax) where T<:Number = Zeros{T}(ax) -mulzeros(::Type{T}, M) where T<:Number = fillzeros(T, axes(M)) -mulzeros(::Type{T}, M::Mul{<:DualLayout,<:Any,<:Adjoint}) where T<:Number = fillzeros(T, axes(M,2))' -mulzeros(::Type{T}, M::Mul{<:DualLayout,<:Any,<:Transpose}) where T<:Number = transpose(fillzeros(T, axes(M,2))) - # initiate array-valued MulAdd function _mulzeros!(dest::AbstractVector{T}, A, B) where T for k in axes(dest,1) From ebab13d7767b713441ad904011eff7dc0818f055 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 4 Aug 2023 21:48:53 +0200 Subject: [PATCH 2/7] import mulzeros --- ext/ArrayLayoutsFillArraysExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/ArrayLayoutsFillArraysExt.jl b/ext/ArrayLayoutsFillArraysExt.jl index 14e16be..ff664fa 100644 --- a/ext/ArrayLayoutsFillArraysExt.jl +++ b/ext/ArrayLayoutsFillArraysExt.jl @@ -5,8 +5,8 @@ using FillArrays: AbstractFill, getindex_value using ArrayLayouts using ArrayLayouts: OnesLayout, Mul, MulAdd -import ArrayLayouts: MemoryLayout, _copyto!, sub_materialize, diagonaldata -export layoutfillmul, mulzeros +import ArrayLayouts: MemoryLayout, _copyto!, sub_materialize, diagonaldata, mulzeros +export layoutfillmul import Base: copy, * import Base.Broadcast: materialize! From cf5407e1a6f917753bd582e7dfa967fb9bafb7ed Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 4 Aug 2023 21:56:52 +0200 Subject: [PATCH 3/7] use diagonal --- ext/ArrayLayoutsFillArraysExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/ArrayLayoutsFillArraysExt.jl b/ext/ArrayLayoutsFillArraysExt.jl index ff664fa..4dff5da 100644 --- a/ext/ArrayLayoutsFillArraysExt.jl +++ b/ext/ArrayLayoutsFillArraysExt.jl @@ -4,7 +4,7 @@ using FillArrays using FillArrays: AbstractFill, getindex_value using ArrayLayouts -using ArrayLayouts: OnesLayout, Mul, MulAdd +using ArrayLayouts: OnesLayout, Mul, MulAdd, diagonal import ArrayLayouts: MemoryLayout, _copyto!, sub_materialize, diagonaldata, mulzeros export layoutfillmul From 34bd574be790615f50758230887a6a32ef1d8442 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 4 Aug 2023 22:11:57 +0200 Subject: [PATCH 4/7] move uniformscaling to extension --- ext/ArrayLayoutsFillArraysExt.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/ext/ArrayLayoutsFillArraysExt.jl b/ext/ArrayLayoutsFillArraysExt.jl index 4dff5da..2e70737 100644 --- a/ext/ArrayLayoutsFillArraysExt.jl +++ b/ext/ArrayLayoutsFillArraysExt.jl @@ -8,7 +8,7 @@ using ArrayLayouts: OnesLayout, Mul, MulAdd, diagonal import ArrayLayouts: MemoryLayout, _copyto!, sub_materialize, diagonaldata, mulzeros export layoutfillmul -import Base: copy, * +import Base: copy, *, +, - import Base.Broadcast: materialize! import LinearAlgebra using LinearAlgebra: Adjoint, Transpose, Symmetric, Hermitian, Diagonal, @@ -92,6 +92,16 @@ sub_materialize(::OnesLayout, V, ax) = Ones{eltype(V)}(ax) _apply(_, _, op, Λ::UniformScaling, A::AbstractMatrix) = op(Diagonal(Fill(Λ.λ,(axes(A,1),))), A) _apply(_, _, op, A::AbstractMatrix, Λ::UniformScaling) = op(A, Diagonal(Fill(Λ.λ,(axes(A,1),)))) +for Typ in (:LayoutMatrix, :(Symmetric{<:Any,<:LayoutMatrix}), :(Hermitian{<:Any,<:LayoutMatrix}), + :(Adjoint{<:Any,<:LayoutMatrix}), :(Transpose{<:Any,<:LayoutMatrix})) + @eval begin + @inline +(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), +, A, Λ) + @inline +(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), +, Λ, A) + @inline -(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), -, A, Λ) + @inline -(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), -, Λ, A) + end +end + # equivalent to rescaling function materialize!(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout}}) checkdimensions(M) From 86e5aae3933f7c46fa4251b7a5ac1cacffe8923f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 4 Aug 2023 22:20:11 +0200 Subject: [PATCH 5/7] remove left-overs --- src/mul.jl | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/mul.jl b/src/mul.jl index 6c1f749..812fb9d 100644 --- a/src/mul.jl +++ b/src/mul.jl @@ -364,11 +364,6 @@ LinearAlgebra.dot(x::AbstractVector, A::Symmetric{<:Real,<:LayoutMatrix}, y::Abs for Typ in (:LayoutMatrix, :(Symmetric{<:Any,<:LayoutMatrix}), :(Hermitian{<:Any,<:LayoutMatrix}), :(Adjoint{<:Any,<:LayoutMatrix}), :(Transpose{<:Any,<:LayoutMatrix})) - @eval begin - @inline Base.power_by_squaring(A::$Typ, p::Integer) = _power_by_squaring(MemoryLayout(A), size(A), A, p) - @inline +(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), +, A, Λ) - @inline +(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), +, Λ, A) - @inline -(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), -, A, Λ) - @inline -(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), -, Λ, A) - end + @eval @inline Base.power_by_squaring(A::$Typ, p::Integer) = + _power_by_squaring(MemoryLayout(A), size(A), A, p) end From 227b95404fa30251743d4544c0d369774b168d57 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 4 Aug 2023 22:28:35 +0200 Subject: [PATCH 6/7] rm uniformscaling addition/subtraction specializations? --- ext/ArrayLayoutsFillArraysExt.jl | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/ext/ArrayLayoutsFillArraysExt.jl b/ext/ArrayLayoutsFillArraysExt.jl index 2e70737..fda7b16 100644 --- a/ext/ArrayLayoutsFillArraysExt.jl +++ b/ext/ArrayLayoutsFillArraysExt.jl @@ -88,20 +88,6 @@ sub_materialize(::OnesLayout, V, ax) = Ones{eltype(V)}(ax) @inline LinearAlgebra.dot(a::LayoutVector, b::AbstractFill{<:Any,1}) = FillArrays._fill_dot_rev(a,b) @inline LinearAlgebra.dot(a::AbstractFill{<:Any,1}, b::LayoutVector) = FillArrays._fill_dot(a,b) -# TODO: Remove unnecessary _apply -_apply(_, _, op, Λ::UniformScaling, A::AbstractMatrix) = op(Diagonal(Fill(Λ.λ,(axes(A,1),))), A) -_apply(_, _, op, A::AbstractMatrix, Λ::UniformScaling) = op(A, Diagonal(Fill(Λ.λ,(axes(A,1),)))) - -for Typ in (:LayoutMatrix, :(Symmetric{<:Any,<:LayoutMatrix}), :(Hermitian{<:Any,<:LayoutMatrix}), - :(Adjoint{<:Any,<:LayoutMatrix}), :(Transpose{<:Any,<:LayoutMatrix})) - @eval begin - @inline +(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), +, A, Λ) - @inline +(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), +, Λ, A) - @inline -(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), -, A, Λ) - @inline -(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), -, Λ, A) - end -end - # equivalent to rescaling function materialize!(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout}}) checkdimensions(M) From 7bbaeb9f91273f4dd8c1146280608d2f7ee329d3 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 23 Aug 2023 11:52:51 +0200 Subject: [PATCH 7/7] Update ArrayLayoutsFillArraysExt.jl --- ext/ArrayLayoutsFillArraysExt.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/ext/ArrayLayoutsFillArraysExt.jl b/ext/ArrayLayoutsFillArraysExt.jl index fda7b16..b8bcc61 100644 --- a/ext/ArrayLayoutsFillArraysExt.jl +++ b/ext/ArrayLayoutsFillArraysExt.jl @@ -39,6 +39,14 @@ end *(a::Adjoint{T, <:LayoutMatrix{T}} where T, b::Zeros{<:Any, 2}) = FillArrays.mult_zeros(a, b) *(A::Adjoint{<:Any, <:Zeros{<:Any,1}}, B::Diagonal{<:Any,<:LayoutVector}) = (B' * A')' *(A::Transpose{<:Any, <:Zeros{<:Any,1}}, B::Diagonal{<:Any,<:LayoutVector}) = transpose(transpose(B) * transpose(A)) +*(a::Adjoint{<:Number,<:LayoutVector}, b::Zeros{<:Number,1})= FillArrays._adjvec_mul_zeros(a, b) +function *(a::Transpose{T, <:LayoutVector{T}}, b::Zeros{T, 1}) where T<:Real + la, lb = length(a), length(b) + if la ≠ lb + throw(DimensionMismatch("dot product arguments have lengths $la and $lb")) + end + return zero(T) +end # equivalent to rescaling function materialize!(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}})