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 by Sebastian Ament <[email protected]>.
renespoerer · Dec 15, 2023 · cb3ecdd · cb3ecdd
1 parent 065ed7e
commit cb3ecdd
Show file tree

Hide file tree

Showing 3 changed files with 53 additions and 21 deletions.
diff --git a/src/ToeplitzMatrices.jl b/src/ToeplitzMatrices.jl
@@ -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")

diff --git a/src/linearalgebra.jl b/src/linearalgebra.jl
@@ -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)
@@ -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])
@@ -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)
@@ -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)
@@ -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
diff --git a/src/toeplitz.jl b/src/toeplitz.jl
@@ -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