Merge 64238d1 into 7ed5934

.github/workflows/CI.yml

@@ -14,7 +14,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1'
           - 'nightly'
         os:

Project.toml

@@ -1,13 +1,23 @@
 name = "ToeplitzMatrices"
 uuid = "c751599d-da0a-543b-9d20-d0a503d91d24"
-version = "0.7.0"
+version = "0.7.1"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 AbstractFFTs = "0.4, 0.5, 1"
+DSP = "0.7"
 StatsBase = "0.32, 0.33"
-julia = "1"
+julia = "1.6"
+
+[extras]
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["FFTW", "Pkg", "Test"]

src/ToeplitzMatrices.jl

@@ -7,12 +7,12 @@ import Base: convert, *, \, getindex, print_matrix, size, Matrix, +, -, copy, si
 import LinearAlgebra: cholesky, cholesky!, eigvals, inv, ldiv!, mul!, pinv, tril, triu
 
 using LinearAlgebra: LinearAlgebra, Adjoint, Factorization, factorize, Cholesky,
-    DimensionMismatch, rmul!
+    DimensionMismatch, rmul!, sym_uplo
 
 using AbstractFFTs
 using AbstractFFTs: Plan
 
-flipdim(A, d) = reverse(A, dims=d)
+using DSP: conv
 
 export Toeplitz, SymmetricToeplitz, Circulant, TriangularToeplitz, Hankel, chan, strang
 
@@ -34,12 +34,9 @@ struct ToeplitzFactorization{T<:Number,A<:AbstractToeplitz{T},S<:Number,P<:Plan{
 end
 
 size(A::AbstractToeplitz) = (size(A, 1), size(A, 2))
-function getindex(A::AbstractToeplitz, i::Integer)
-    return A[mod(i - 1, size(A, 1)) + 1, div(i - 1, size(A, 1)) + 1]
-end
 
-convert(::Type{AbstractMatrix{T}}, S::AbstractToeplitz) where {T} = convert(AbstractToeplitz{T}, S)
-convert(::Type{AbstractArray{T}}, S::AbstractToeplitz) where {T} = convert(AbstractToeplitz{T}, S)
+convert(::Type{AbstractMatrix{T}}, S::AbstractToeplitz) where {T<:Number} = convert(AbstractToeplitz{T}, S)
+convert(::Type{AbstractArray{T}}, S::AbstractToeplitz) where {T<:Number} = convert(AbstractToeplitz{T}, S)
 
 # Fast application of a general Toeplitz matrix to a column vector via FFT
 function mul!(
@@ -248,11 +245,11 @@ function getindex(A::Toeplitz, i::Integer, j::Integer)
 end
 
 # Form a lower triangular Toeplitz matrix by annihilating all entries above the k-th diaganal
-function tril(A::Toeplitz, k = 0)
+function tril(A::Toeplitz, k::Integer = 0)
     if k > 0
         error("Second argument cannot be positive")
     end
-    Al = TriangularToeplitz(copy(A.vc), 'L', length(A.vr))
+    Al = TriangularToeplitz(copy(A.vc), :L)
     if k < 0
         for i in 1:-k
             Al.ve[i] = zero(eltype(A))
@@ -262,11 +259,11 @@ function tril(A::Toeplitz, k = 0)
 end
 
 # Form a lower triangular Toeplitz matrix by annihilating all entries below the k-th diagonal
-function triu(A::Toeplitz, k = 0)
+function triu(A::Toeplitz, k::Integer = 0)
     if k < 0
         error("Second argument cannot be negative")
     end
-    Al = TriangularToeplitz(copy(A.vr), 'U', length(A.vc))
+    Al = TriangularToeplitz(copy(A.vr), :U)
     if k > 0
         for i in 1:k
             Al.ve[i] = zero(eltype(A))
@@ -384,7 +381,8 @@ function getindex(C::Circulant, i::Integer, j::Integer)
     @boundscheck if i < 1 || i > n || j < 1 || j > n
         throw(BoundsError(C, (i,j)))
     end
-    return C.vc[mod(i - j, length(C.vc)) + 1]
+    d = i - j
+    return C.vc[d < 0 ? n+d+1 : d+1]
 end
 
 LinearAlgebra.ldiv!(C::Circulant, b::AbstractVector) = ldiv!(factorize(C), b)
@@ -415,9 +413,9 @@ function Base.inv(C::Circulant)
     vc = F.dft \ vdft
     return Circulant(maybereal(eltype(C), vc))
 end
-function Base.inv(C::CirculantFactorization)
+function Base.inv(C::CirculantFactorization{T,S,P}) where {T,S,P}
     vdft = map(inv, C.vcvr_dft)
-    return CirculantFactorization(vdft, similar(vdft), C.dft)
+    return CirculantFactorization{T,S,P}(vdft, similar(vdft), C.dft)
 end
 
 function strang(A::AbstractMatrix{T}) where T
@@ -490,26 +488,16 @@ function Base.:*(A::CirculantFactorization, B::CirculantFactorization)
     return Circulant(maybereal(Base.promote_eltype(A, B), vc))
 end
 
-Base.:*(A::Adjoint{<:Circulant}, B::Circulant) = factorize(parent(A))' * factorize(B)
-Base.:*(A::Adjoint{<:Circulant}, B::CirculantFactorization) = factorize(parent(A))' * B
-Base.:*(A::Adjoint{<:CirculantFactorization}, B::Circulant) = A * factorize(B)
-function Base.:*(A::Adjoint{<:CirculantFactorization}, B::CirculantFactorization)
-    C = parent(A)
-    C_vcvr_dft = C.vcvr_dft
-    B_vcvr_dft = B.vcvr_dft
-    m = length(C_vcvr_dft)
-    n = length(B_vcvr_dft)
-    if m != n
-        throw(DimensionMismatch(
-            "size of matrix A, $(m)x$(m), does not match size of matrix B, $(n)x$(n)"
-        ))
-    end
-
-    tmp = conj.(C_vcvr_dft) .* B_vcvr_dft
-    vc = C.dft \ tmp
-
-    return Circulant(maybereal(Base.promote_eltype(C, B), vc))
-end
+# Make an eager adjoint, similar to adjoints of Diagonal in LinearAlgebra
+adjoint(C::ToeplitzFactorization{T,Circulant{T},S,P}) where {T,S,P} =
+    ToeplitzFactorization{T,Circulant{T},S,P}(conj.(C.vcvr_dft), C.tmp, C.dft)
+Base.:*(A::Adjoint{<:Any,<:Circulant}, B::Circulant) = factorize(parent(A))' * factorize(B)
+Base.:*(A::Adjoint{<:Any,<:Circulant}, B::CirculantFactorization) =
+    factorize(parent(A))' * B
+Base.:*(A::Adjoint{<:Any,<:Circulant}, B::Adjoint{<:Any,<:Circulant}) =
+    factorize(parent(A))' * factorize(parent(B))'
+Base.:*(A::Circulant, B::Adjoint{<:Any,<:Circulant}) = factorize(A) * factorize(parent(B))'
+Base.:*(A::CirculantFactorization, B::Adjoint{<:Any,<:Circulant}) = A * factorize(parent(B))'
 
 (*)(scalar::Number, C::Circulant) = Circulant(scalar * C.vc)
 (*)(C::Circulant,scalar::Number) = Circulant(C.vc * scalar)
@@ -564,7 +552,7 @@ end
 
 convert(::Type{AbstractToeplitz{T}}, A::TriangularToeplitz) where {T} = convert(TriangularToeplitz{T},A)
 convert(::Type{TriangularToeplitz{T}}, A::TriangularToeplitz) where {T} =
-    TriangularToeplitz(convert(Vector{T},A.ve),A.uplo=='U' ? (:U) : (:L))
+    TriangularToeplitz(convert(Vector{T}, A.ve), A.uplo)
 
 
 function size(A::TriangularToeplitz, dim::Int)
@@ -587,17 +575,17 @@ function getindex(A::TriangularToeplitz{T}, i::Integer, j::Integer) where T
 end
 
 function (*)(A::TriangularToeplitz, B::TriangularToeplitz)
-    n = size(A, 1)
+    n = size(A, 2)
     if n != size(B, 1)
         throw(DimensionMismatch(""))
     end
     if A.uplo == B.uplo
         return TriangularToeplitz(conv(A.ve, B.ve)[1:n], A.uplo)
     end
-    return Triangular(Matrix(A), A.uplo) * Triangular(Matrix(B), B.uplo)
+    return A * Matrix(B)
 end
 
-function Base.:*(A::Adjoint{<:TriangularToeplitz}, b::AbstractVector)
+function Base.:*(A::Adjoint{<:Any,<:TriangularToeplitz}, b::AbstractVector)
     M = parent(A)
     return TriangularToeplitz{eltype(M)}(M.ve, M.uplo) * b
 end
@@ -614,23 +602,27 @@ function smallinv(A::TriangularToeplitz{T}) where T
         end
         b[k] = -tmp/A.ve[1]
     end
-    return TriangularToeplitz(b, symbol(A.uplo))
+    return TriangularToeplitz(b, A.uplo)
 end
 
 function inv(A::TriangularToeplitz{T}) where T
     n = size(A, 1)
     if n <= 64
         return smallinv(A)
     end
-    np2 = nextpow2(n)
+    np2 = nextpow(2, n)
     if n != np2
-        return TriangularToeplitz(inv(TriangularToeplitz([A.ve, zeros(T, np2 - n)],
-            symbol(A.uplo))).ve[1:n], symbol(A.uplo))
+        return TriangularToeplitz(
+            inv(TriangularToeplitz([A.ve; zeros(T, np2 - n)], sym_uplo(A.uplo))).ve[1:n],
+            sym_uplo(A.uplo)
+        )
     end
     nd2 = div(n, 2)
-    a1 = inv(TriangularToeplitz(A.ve[1:nd2], symbol(A.uplo))).ve
-    return TriangularToeplitz([a1, -(TriangularToeplitz(a1, symbol(A.uplo)) *
-        (Toeplitz(A.ve[nd2 + 1:end], A.ve[nd2 + 1:-1:2]) * a1))], symbol(A.uplo))
+    a1 = inv(TriangularToeplitz(A.ve[1:nd2], sym_uplo(A.uplo))).ve
+    return TriangularToeplitz(
+        [a1; -(TriangularToeplitz(a1, sym_uplo(A.uplo)) * (Toeplitz(A.ve[nd2 + 1:end], A.ve[nd2 + 1:-1:2]) * a1))],
+        sym_uplo(A.uplo)
+    )
 end
 
 # ldiv!(A::TriangularToeplitz,b::StridedVector) = inv(A)*b
@@ -746,9 +738,13 @@ end
 
 # Fast application of a general Hankel matrix to a general vector
 *(A::Hankel, b::AbstractVector) = A.T * reverse(b)
-
+mul!(y::StridedVector, A::Hankel, x::StridedVector, α::Number, β::Number) =
+    mul!(y, A.T, view(x, reverse(axes(x, 1))), α, β)
 # Fast application of a general Hankel matrix to a general matrix
-*(A::Hankel, B::AbstractMatrix) = A.T * flipdim(B, 1)
+*(A::Hankel, B::AbstractMatrix) = A.T * reverse(B, dims=1)
+mul!(Y::StridedMatrix, A::Hankel, X::StridedMatrix, α::Number, β::Number) =
+    mul!(Y, A.T, view(X, reverse(axes(X, 1)), :), α, β)
+
 ## BigFloat support
 
 (*)(A::Toeplitz{T}, b::AbstractVector) where {T<:BigFloat} = irfft(

test/runtests.jl

@@ -9,15 +9,14 @@ end
 
 using ToeplitzMatrices, StatsBase, Test, LinearAlgebra
 
-using Base: copyto!
 using FFTW: fft
 
 ns = 101
 nl = 2000
 
 xs = randn(ns, 5)
 xl = randn(nl, 5)
-vc = LinRange(1,3,3) # for testing with AbstractVector
+vc = 1.0:3.0 # for testing with AbstractVector
 vv = Vector(vc)
 vr = [1, 5.]
 
@@ -119,13 +118,13 @@ end
         @test diag(H) == [1,3,5,7,9]
 
         x = ones(5)
-        @test Matrix(H)*x ≈ H*x
+        @test mul!(copy(x), H, x) ≈ Matrix(H)*x ≈ H*x
 
         Hs = Hankel(0.9.^(ns-1:-1:0), 0.4.^(0:ns-1))
         Hl = Hankel(0.9.^(nl-1:-1:0), 0.4.^(0:nl-1))
         @test Hs * xs[:,1] ≈ Matrix(Hs) * xs[:,1]
-        @test Hs * xs ≈ Matrix(Hs) * xs
-        @test Hl * xl ≈ Matrix(Hl) * xl
+        @test mul!(copy(xs), Hs, xs) ≈ Hs * xs ≈ Matrix(Hs) * xs
+        @test mul!(copy(xl), Hl, xl) ≈ Hl * xl ≈ Matrix(Hl) * xl
         @test Matrix(Hankel(reverse(vc),vr)) == Matrix(Hankel(reverse(vv),vr))
     end
 
@@ -203,6 +202,7 @@ end
     @test isa(convert(AbstractArray{ComplexF64},T),TriangularToeplitz{ComplexF64})
     @test isa(convert(ToeplitzMatrices.AbstractToeplitz{ComplexF64},T),TriangularToeplitz{ComplexF64})
     @test isa(convert(ToeplitzMatrices.TriangularToeplitz{ComplexF64},T),TriangularToeplitz{ComplexF64})
+    @test isa(convert(Toeplitz, T), Toeplitz)
 
     T = Hankel(ones(2),ones(2))
 
@@ -259,9 +259,11 @@ end
     M4 = Matrix(C4)
     M5 = Matrix(C5)
 
-    C = C1*C2
-    @test C isa Circulant
-    @test C ≈ M1*M2
+    for t1 in (identity, adjoint), t2 in (identity, adjoint), fact in (identity, factorize)
+        C = t1(fact(C1))*t2(fact(C2))
+        @test C isa Circulant
+        @test C ≈ t1(M1)*t2(M2)
+    end
 
     C = C1-C2
     @test C isa Circulant
@@ -330,9 +332,39 @@ end
 
     # Test for issue #47
     I = inv(C1)*C1
+    I2 = inv(factorize(C1))*C1
     e = rand(5)
     # I should be close to identity
-    @test I*e ≈ e
+    @test I*e ≈ I2*e ≈ e
+end
+
+@testset "TriangularToeplitz" begin
+    A = [1.0 2.0;
+         3.0 4.0]
+    TU = TriangularToeplitz(A, :U)
+    TL = TriangularToeplitz(A, :L)
+    @test (TU * TU)::TriangularToeplitz ≈ Matrix(TU)*Matrix(TU)
+    @test (TL * TL)::TriangularToeplitz ≈ Matrix(TL)*Matrix(TL)
+    @test (TU * TL) ≈ Matrix(TU)*Matrix(TL)
+    for T in (TU, TL)
+        @test inv(T)::TriangularToeplitz ≈ inv(Matrix(T))
+    end
+    T = Toeplitz(A)
+    TU = triu(T)
+    @test TU isa TriangularToeplitz
+    @test istriu(TU)
+    @test TU == Toeplitz(triu(A))
+    TL = tril(T)
+    @test TL isa TriangularToeplitz
+    @test istril(TL)
+    @test TL == Toeplitz(tril(A))
+    for n in (65, 128)
+        A = randn(n, n)
+        TU = TriangularToeplitz(A, :U)
+        TL = TriangularToeplitz(A, :L)
+        @test_broken inv(TU)::TriangularToeplitz ≈ inv(Matrix(TU))
+        @test inv(TL)::TriangularToeplitz ≈ inv(Matrix(TL))
+    end
 end
 
 @testset "Cholesky" begin