Skip to content

Commit

Permalink
added base functions for ToeplitzFactorization
Browse files Browse the repository at this point in the history
Added adjoint, copy, size, and length for ToeplitzFactorization; added
fields n and m to facilitate this.

Copied from PR JuliaLinearAlgebra/pull/65 by
Sebastian Ament <sebastianeament@gmail.com>.
  • Loading branch information
renespoerer committed Dec 15, 2023
1 parent 065ed7e commit cb3ecdd
Show file tree
Hide file tree
Showing 3 changed files with 53 additions and 21 deletions.
11 changes: 0 additions & 11 deletions src/ToeplitzMatrices.jl
Original file line number Diff line number Diff line change
Expand Up @@ -72,17 +72,6 @@ function isdiag(A::AbstractToeplitz)
all(iszero, @view vr[2:end]) && all(iszero, @view vc[2:end])
end

"""
ToeplitzFactorization
Factorization of a Toeplitz matrix using FFT.
"""
struct ToeplitzFactorization{T,A<:AbstractToeplitz{T},S<:Number,P<:Plan{S}} <: Factorization{T}
vcvr_dft::Vector{S}
tmp::Vector{S}
dft::P
end

include("toeplitz.jl")
include("special.jl")
include("hankel.jl")
Expand Down
21 changes: 11 additions & 10 deletions src/linearalgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -121,13 +121,13 @@ end

function factorize(A::Toeplitz)
T = eltype(A)
m, n = size(A)
S = promote_type(float(T), Complex{Float32})
tmp = Vector{S}(undef, m + n - 1)
n, m = size(A)
S = promote_type(T, Complex{Float32})
tmp = Vector{S}(undef, n + m - 1)
copyto!(tmp, A.vc)
copyto!(tmp, m + 1, Iterators.reverse(A.vr), 1, n - 1)
copyto!(tmp, n + 1, Iterators.reverse(A.vr), 1, m - 1)
dft = plan_fft!(tmp)
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft)
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, n, m)
end

function ldiv!(A::Toeplitz, b::StridedVector)
Expand All @@ -146,7 +146,7 @@ function factorize(A::SymmetricToeplitz{T}) where {T<:Number}
@inbounds tmp[m + 1] = zero(T)
copyto!(tmp, m + 2, Iterators.reverse(vc), 1, m - 1)
dft = plan_fft!(tmp)
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft)
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, m, m)
end

ldiv!(A::SymmetricToeplitz, b::StridedVector) = copyto!(b, IterativeLinearSolvers.cg(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1])
Expand Down Expand Up @@ -184,11 +184,12 @@ const CirculantFactorization{T, V<:AbstractVector{T}} = ToeplitzFactorization{T,
function factorize(C::Circulant)
T = eltype(C)
vc = C.vc
m = length(vc)
S = promote_type(float(T), Complex{Float32})
tmp = Vector{S}(undef, length(vc))
tmp = Vector{S}(undef, m)
copyto!(tmp, vc)
dft = plan_fft!(tmp)
return ToeplitzFactorization{T,typeof(C),S,typeof(dft)}(dft * tmp, similar(tmp), dft)
return ToeplitzFactorization{T,typeof(C),S,typeof(dft)}(dft * tmp, similar(tmp), dft, m, m)
end

Base.:*(A::Circulant, B::Circulant) = factorize(A) * factorize(B)
Expand Down Expand Up @@ -420,7 +421,7 @@ function factorize(A::LowerTriangularToeplitz)
tmp = zeros(S, 2 * n - 1)
copyto!(tmp, v)
dft = plan_fft!(tmp)
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft)
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, n, n)
end
function factorize(A::UpperTriangularToeplitz)
T = eltype(A)
Expand All @@ -431,5 +432,5 @@ function factorize(A::UpperTriangularToeplitz)
tmp[1] = v[1]
copyto!(tmp, n + 1, Iterators.reverse(v), 1, n - 1)
dft = plan_fft!(tmp)
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft)
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, n, n)
end
42 changes: 42 additions & 0 deletions src/toeplitz.jl
Original file line number Diff line number Diff line change
Expand Up @@ -144,3 +144,45 @@ function rmul!(A::Toeplitz, x::Number)
rmul!(A.vr, x)
A
end

"""
ToeplitzFactorization
Factorization of a Toeplitz matrix using FFT.
"""
struct ToeplitzFactorization{T,A<:AbstractToeplitz{T},S<:Number,P<:Plan{S}} <: Factorization{T}
vcvr_dft::Vector{S}
tmp::Vector{S}
dft::P
n::Int
m::Int
end

# doing this non-lazily simplifies implementation of mul!, ldiv! for adjoints
# of Toeplitz factorizations significantly Base.adjoint(A::ToeplitzFactorization) = Adjoint(A)
adjoint(T::ToeplitzFactorization) = adjoint!(copy(T))
function copy(T::ToeplitzFactorization)
vcvr_dft = copy(T.vcvr_dft)
dft = plan_fft!(vcvr_dft)
typeof(T)(vcvr_dft, copy(T.tmp), dft, T.n, T.m)
end
# calculates the adjoint but reuses the temporary memory of T
function adjoint!(T::ToeplitzFactorization)
@. T.vcvr_dft = conj(T.vcvr_dft)
typeof(T)(T.vcvr_dft, T.tmp, T.dft, T.m, T.n) # switching n and m
end

size(A::ToeplitzFactorization) = (A.n, A.m)
function Base.size(A::ToeplitzFactorization, i::Int)
if i == 1
A.n
elseif i == 2
A.m
elseif i > 2
1
else
throw(DomainError("dimension i cannot be non-positive"))
end
end

length(A::ToeplitzFactorization) = A.n * A.m

0 comments on commit cb3ecdd

Please sign in to comment.