diff --git a/README.md b/README.md index b5da0b7..f9847f0 100644 --- a/README.md +++ b/README.md @@ -82,10 +82,7 @@ The `reverse` can transform between Hankel and Toeplitz. It is used to achieve f |copy|✓|✓|✓|✓|✓|✓| |similar|✓|✓|✓|✓|✓|✓| |zero|✓|✓|✓|✓|✓|✓| -|real|✓|✓|✓|✓|✓|✓| -|imag|✓|✓|✓|✓|✓|✓| |fill!|✓|✗|✗|✗|✗|✓| -|conj|✓|✓|✓|✓|✓|✓| |transpose|✓|✓|✓|✓|✓|✓| |adjoint|✓|✓|✓|✓|✓|✓| |tril!|✓|✗|✗|✓|✓|✗| @@ -94,22 +91,22 @@ The `reverse` can transform between Hankel and Toeplitz. It is used to achieve f |triu|✓|✓|✓|✓|✓|✗| |+|✓|✓|✓|✓|✓|✓| |-|✓|✓|✓|✓|✓|✓| -|scalar
mult|✓|✓|✓|✓|✓|✓| |==|✓|✓|✓|✓|✓|✓| |issymmetric||||||| |istriu||||||| |istril||||||| -|iszero|✓|✓|✓|✓|✓|| |isone||||||| |diag|✓|✓|✓|✓|✓|✓| |copyto!|✓|✓|✓|✓|✓|✓| |reverse|✓|✓|✓|✓|✓|✓| -|broadcast||||||| +|\_all|✓|✓|✓|✓|✓|✓| +|unary broadcast|✓|✓|✓|✓|✓|✓| +|number broadcast|✓|✓|✓|✓|✓|✓| |broadcast!||||||| -Note that scalar multiplication, `conj`, `+` and `-` could be removed once `broadcast` is implemented. +Note that `+` and `-` could be removed once binary `broadcast` is implemented. `reverse(Hankel)` returns a `Toeplitz`, while `reverse(AbstractToeplitz)` returns a `Hankel`. diff --git a/src/ToeplitzMatrices.jl b/src/ToeplitzMatrices.jl index 5edcd81..fd5730b 100644 --- a/src/ToeplitzMatrices.jl +++ b/src/ToeplitzMatrices.jl @@ -1,10 +1,12 @@ module ToeplitzMatrices import DSP: conv -import Base: adjoint, convert, transpose, size, getindex, similar, copy, getproperty, inv, sqrt, copyto!, reverse, conj, zero, fill!, checkbounds, real, imag, isfinite, DimsInteger, iszero +import Base: convert, size, getindex, similar, copy, getproperty, copyto!, reverse, zero, fill!, checkbounds, isfinite, DimsInteger, iszero +import Base: all, adjoint, transpose, real, imag, inv, sqrt, conj import Base: parent import Base: ==, +, -, *, \ import Base: AbstractMatrix +import Base.Broadcast: broadcasted, DefaultMatrixStyle import LinearAlgebra: Cholesky, Factorization import LinearAlgebra: ldiv!, factorize, lmul!, pinv, eigvals, eigvecs, eigen, Eigen, det import LinearAlgebra: cholesky!, cholesky, tril!, triu!, checksquare, rmul!, dot, mul!, tril, triu, diag @@ -34,7 +36,7 @@ end include("iterativeLinearSolvers.jl") # Abstract -abstract type AbstractToeplitz{T<:Number} <: AbstractMatrix{T} end +abstract type AbstractToeplitz{T} <: AbstractMatrix{T} end size(A::AbstractToeplitz) = (length(A.vc),length(A.vr)) @inline _vr(A::AbstractToeplitz) = A.vr @@ -46,6 +48,11 @@ AbstractArray{T}(A::AbstractToeplitz) where T = AbstractToeplitz{T}(A) convert(::Type{AbstractToeplitz{T}}, A::AbstractToeplitz) where T = AbstractToeplitz{T}(A) isconcrete(A::AbstractToeplitz) = isconcretetype(typeof(A.vc)) && isconcretetype(typeof(A.vr)) + +for op in (:+, :max) + @eval mapreduce(f, ::typeof($op), A::AbstractToeplitz) = mapreduce(identity, $op, [mapreduce(f, $op, diag(A, k)) for k in -size(A,1)+1:size(A,2)-1]) +end +all(f, A::AbstractToeplitz, ::Colon) = all(f, A.vc, :) && all(f, A.vr, :) iszero(A::AbstractToeplitz) = iszero(A.vc) && iszero(A.vr) function diag(A::AbstractToeplitz, k::Integer=0) if k >= size(A, 2) || -k >= size(A, 1) diff --git a/src/hankel.jl b/src/hankel.jl index c3e396b..ab25cc7 100644 --- a/src/hankel.jl +++ b/src/hankel.jl @@ -60,6 +60,10 @@ end convert(::Type{AbstractArray{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A) convert(::Type{AbstractMatrix{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A) convert(::Type{Hankel{T}}, A::Hankel) where {T<:Number} = Hankel{T}(convert(AbstractVector{T}, A.v), size(A)) +broadcasted(::DefaultMatrixStyle, f, A::Hankel) = Hankel(f.(A.v), A.size) +broadcasted(::DefaultMatrixStyle, f, x::Number, A::Hankel) = Hankel(f.(x, A.v), A.size) +broadcasted(::DefaultMatrixStyle, f, A::Hankel, x::Number) = Hankel(f.(A.v, x), A.size) +all(f, A::Hankel, ::Colon) = all(f, A.v, :) # Size size(H::Hankel) = H.size @@ -72,7 +76,7 @@ Base.@propagate_inbounds function getindex(A::Hankel, i::Integer, j::Integer) return A.v[i+j-1] end AbstractMatrix{T}(A::Hankel) where {T} = Hankel{T}(AbstractVector{T}(A.v), A.size) -for fun in (:zero, :conj, :copy, :-, :similar, :real, :imag) +for fun in (:zero, :copy) @eval $fun(A::Hankel) = Hankel($fun(A.v), size(A)) end for op in (:+, :-) @@ -102,8 +106,6 @@ end transpose(A::Hankel) = Hankel(A.v, reverse(size(A))) adjoint(A::Hankel) = transpose(conj(A)) (==)(A::Hankel, B::Hankel) = A.v == B.v && size(A) == size(B) -(*)(scalar::Number, C::Hankel) = Hankel(scalar * C.v, size(C)) -(*)(C::Hankel,scalar::Number) = Hankel(C.v * scalar, size(C)) isconcrete(A::Hankel) = isconcretetype(A.v) diff --git a/src/special.jl b/src/special.jl index f01d790..9f1675d 100644 --- a/src/special.jl +++ b/src/special.jl @@ -74,7 +74,8 @@ for TYPE in (:SymmetricToeplitz, :Circulant, :LowerTriangularToeplitz, :UpperTri @eval $TYPE(v::$TY) = $TYPE{eltype(v)}(v) end end -TriangularToeplitz{T}=Union{UpperTriangularToeplitz{T},LowerTriangularToeplitz{T}} +SymToeplitzOrCirc{T} = Union{SymmetricToeplitz{T}, Circulant{T}} +TriangularToeplitz{T} = Union{UpperTriangularToeplitz{T}, LowerTriangularToeplitz{T}} # vc and vr function getproperty(A::SymmetricToeplitz, s::Symbol) @@ -129,6 +130,22 @@ transpose(A::Circulant) = Circulant(A.vr) transpose(A::LowerTriangularToeplitz) = UpperTriangularToeplitz(A.v) transpose(A::UpperTriangularToeplitz) = LowerTriangularToeplitz(A.v) +# _all +all(f, A::SymToeplitzOrCirc, ::Colon) = all(f, A.v, :) +all(f, A::TriangularToeplitz, ::Colon) = f(zero(eltype(A))) && all(f, A.v, :) + +# broadcast +for TYPE in (:SymmetricToeplitz, :Circulant) + @eval broadcasted(::DefaultMatrixStyle, f, A::$TYPE) = $TYPE(f.(A.v)) + @eval broadcasted(::DefaultMatrixStyle, f, x::Number, A::$TYPE) = $TYPE(f.(x, A.v)) + @eval broadcasted(::DefaultMatrixStyle, f, A::$TYPE, x::Number) = $TYPE(f.(A.v, x)) +end +for TYPE in (:UpperTriangularToeplitz, :LowerTriangularToeplitz) + @eval broadcasted(::DefaultMatrixStyle, f, A::$TYPE) = iszero(f(zero(eltype(A)))) ? $TYPE(f.(A.v)) : _toep_broadcast(f, A) + @eval broadcasted(::DefaultMatrixStyle, f, x::Number, A::$TYPE) = iszero(f(x, zero(eltype(A)))) ? $TYPE(f.(x, A.v)) : _toep_broadcast(f, x, A) + @eval broadcasted(::DefaultMatrixStyle, f, A::$TYPE, x::Number) = iszero(f(zero(eltype(A)), x)) ? $TYPE(f.(A.v, x)) : _toep_broadcast(f, A, x) +end + # getindex Base.@propagate_inbounds function getindex(A::SymmetricToeplitz, i::Integer, j::Integer) @boundscheck checkbounds(A,i,j) diff --git a/src/toeplitz.jl b/src/toeplitz.jl index 8c98979..ca3216f 100644 --- a/src/toeplitz.jl +++ b/src/toeplitz.jl @@ -10,7 +10,7 @@ struct Toeplitz{T, VC<:AbstractVector{T}, VR<:AbstractVector{T}} <: AbstractToep function Toeplitz{T, VC, VR}(vc::VC, vr::VR) where {T, VC<:AbstractVector{T}, VR<:AbstractVector{T}} require_one_based_indexing(vr, vc) - if first(vc) != first(vr) + if !isequal(first(vc), first(vr)) error("First element of the vectors must be the same") end return new{T,VC,VR}(vc, vr) @@ -54,6 +54,13 @@ Base.@propagate_inbounds function getindex(A::AbstractToeplitz, i::Integer, j::I end end +broadcasted(::DefaultMatrixStyle, f, A::AbstractToeplitz) = _toep_broadcast(f, A) +broadcasted(::DefaultMatrixStyle, f, x::Number, A::AbstractToeplitz) = _toep_broadcast(f, x, A) +broadcasted(::DefaultMatrixStyle, f, A::AbstractToeplitz, x::Number) = _toep_broadcast(f, A, x) +_toep_broadcast(f, A::AbstractToeplitz) = Toeplitz(f.(A.vc), f.(A.vr)) +_toep_broadcast(f, x::Number, A::AbstractToeplitz) = Toeplitz(f.(x, A.vc), f.(x, A.vr)) +_toep_broadcast(f, A::AbstractToeplitz, x::Number) = Toeplitz(f.(A.vc, x), f.(A.vr, x)) + checknonaliased(A::Toeplitz) = Base.mightalias(A.vc, A.vr) && throw(ArgumentError("Cannot modify Toeplitz matrices in place with aliased data")) function tril!(A::Toeplitz, k::Integer=0) @@ -110,7 +117,7 @@ function AbstractMatrix{T}(A::AbstractToeplitz) where {T} vr = AbstractVector{T}(_vr(A)) Toeplitz{T}(vc,vr) end -for fun in (:zero, :conj, :copy, :-, :real, :imag) +for fun in (:zero, :copy) @eval $fun(A::AbstractToeplitz)=Toeplitz($fun(A.vc),$fun(A.vr)) end for op in (:+, :-) @@ -129,8 +136,6 @@ function fill!(A::Toeplitz, x::Number) fill!(A.vr,x) A end -(*)(scalar::Number, C::AbstractToeplitz) = Toeplitz(scalar * C.vc, scalar * C.vr) -(*)(C::AbstractToeplitz,scalar::Number) = Toeplitz(C.vc * scalar, C.vr * scalar) function lmul!(x::Number, A::Toeplitz) checknonaliased(A) diff --git a/test/runtests.jl b/test/runtests.jl index 42eaae3..c1f0c55 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -430,6 +430,23 @@ end end end +@testset "Broadcast" begin + A = rand(ComplexF64, 3, 3) + T = Toeplitz(A) + S = Symmetric(T) + C = Circulant(T) + U = UpperTriangular(T) + L = LowerTriangular(T) + + for M in (T, S, C, U, L) + @testset "$(typeof(M))" begin + @test M .+ 1 == Matrix(M) .+ 1 + @test exp.(M) == exp.(Matrix(M)) + @test isa(M .* 2, typeof(M)) + end + end +end + @testset "Circulant mathematics" begin C1 = Circulant(rand(5)) C2 = Circulant(rand(5))