Skip to content

Commit 276ae92

Browse files
added base functions for ToeplitzFactorization
Added adjoint, copy, size, and length for ToeplitzFactorization; added fields n and m to facilitate this. Copied from PR JuliaLinearAlgebra/pull/65 Co-authored-by: Sebastian Ament <[email protected]>
1 parent 065ed7e commit 276ae92

File tree

3 files changed

+53
-21
lines changed

3 files changed

+53
-21
lines changed

src/ToeplitzMatrices.jl

Lines changed: 0 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -72,17 +72,6 @@ function isdiag(A::AbstractToeplitz)
7272
all(iszero, @view vr[2:end]) && all(iszero, @view vc[2:end])
7373
end
7474

75-
"""
76-
ToeplitzFactorization
77-
78-
Factorization of a Toeplitz matrix using FFT.
79-
"""
80-
struct ToeplitzFactorization{T,A<:AbstractToeplitz{T},S<:Number,P<:Plan{S}} <: Factorization{T}
81-
vcvr_dft::Vector{S}
82-
tmp::Vector{S}
83-
dft::P
84-
end
85-
8675
include("toeplitz.jl")
8776
include("special.jl")
8877
include("hankel.jl")

src/linearalgebra.jl

Lines changed: 11 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -121,13 +121,13 @@ end
121121

122122
function factorize(A::Toeplitz)
123123
T = eltype(A)
124-
m, n = size(A)
125-
S = promote_type(float(T), Complex{Float32})
126-
tmp = Vector{S}(undef, m + n - 1)
124+
n, m = size(A)
125+
S = promote_type(T, Complex{Float32})
126+
tmp = Vector{S}(undef, n + m - 1)
127127
copyto!(tmp, A.vc)
128-
copyto!(tmp, m + 1, Iterators.reverse(A.vr), 1, n - 1)
128+
copyto!(tmp, n + 1, Iterators.reverse(A.vr), 1, m - 1)
129129
dft = plan_fft!(tmp)
130-
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft)
130+
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, n, m)
131131
end
132132

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

152152
ldiv!(A::SymmetricToeplitz, b::StridedVector) = copyto!(b, IterativeLinearSolvers.cg(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1])
@@ -184,11 +184,12 @@ const CirculantFactorization{T, V<:AbstractVector{T}} = ToeplitzFactorization{T,
184184
function factorize(C::Circulant)
185185
T = eltype(C)
186186
vc = C.vc
187+
m = length(vc)
187188
S = promote_type(float(T), Complex{Float32})
188-
tmp = Vector{S}(undef, length(vc))
189+
tmp = Vector{S}(undef, m)
189190
copyto!(tmp, vc)
190191
dft = plan_fft!(tmp)
191-
return ToeplitzFactorization{T,typeof(C),S,typeof(dft)}(dft * tmp, similar(tmp), dft)
192+
return ToeplitzFactorization{T,typeof(C),S,typeof(dft)}(dft * tmp, similar(tmp), dft, m, m)
192193
end
193194

194195
Base.:*(A::Circulant, B::Circulant) = factorize(A) * factorize(B)
@@ -420,7 +421,7 @@ function factorize(A::LowerTriangularToeplitz)
420421
tmp = zeros(S, 2 * n - 1)
421422
copyto!(tmp, v)
422423
dft = plan_fft!(tmp)
423-
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft)
424+
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, n, n)
424425
end
425426
function factorize(A::UpperTriangularToeplitz)
426427
T = eltype(A)
@@ -431,5 +432,5 @@ function factorize(A::UpperTriangularToeplitz)
431432
tmp[1] = v[1]
432433
copyto!(tmp, n + 1, Iterators.reverse(v), 1, n - 1)
433434
dft = plan_fft!(tmp)
434-
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft)
435+
return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, n, n)
435436
end

src/toeplitz.jl

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -144,3 +144,45 @@ function rmul!(A::Toeplitz, x::Number)
144144
rmul!(A.vr, x)
145145
A
146146
end
147+
148+
"""
149+
ToeplitzFactorization
150+
151+
Factorization of a Toeplitz matrix using FFT.
152+
"""
153+
struct ToeplitzFactorization{T,A<:AbstractToeplitz{T},S<:Number,P<:Plan{S}} <: Factorization{T}
154+
vcvr_dft::Vector{S}
155+
tmp::Vector{S}
156+
dft::P
157+
n::Int
158+
m::Int
159+
end
160+
161+
# doing this non-lazily simplifies implementation of mul!, ldiv! for adjoints
162+
# of Toeplitz factorizations significantly Base.adjoint(A::ToeplitzFactorization) = Adjoint(A)
163+
adjoint(T::ToeplitzFactorization) = adjoint!(copy(T))
164+
function copy(T::ToeplitzFactorization)
165+
vcvr_dft = copy(T.vcvr_dft)
166+
dft = plan_fft!(vcvr_dft)
167+
typeof(T)(vcvr_dft, copy(T.tmp), dft, T.n, T.m)
168+
end
169+
# calculates the adjoint but reuses the temporary memory of T
170+
function adjoint!(T::ToeplitzFactorization)
171+
@. T.vcvr_dft = conj(T.vcvr_dft)
172+
typeof(T)(T.vcvr_dft, T.tmp, T.dft, T.m, T.n) # switching n and m
173+
end
174+
175+
size(A::ToeplitzFactorization) = (A.n, A.m)
176+
function Base.size(A::ToeplitzFactorization, i::Int)
177+
if i == 1
178+
A.n
179+
elseif i == 2
180+
A.m
181+
elseif i > 2
182+
1
183+
else
184+
throw(DomainError("dimension i cannot be non-positive"))
185+
end
186+
end
187+
188+
length(A::ToeplitzFactorization) = A.n * A.m

0 commit comments

Comments
 (0)