Unclear how to use pre-computed factorization

[adomasbaliuka]

The README suggests that when using the same matrix multiple times, precomputing the factorization (which uses FFT) can be useful for better performance. However, using this factorization directly does not seem to be supported:

julia> using ToeplitzMatrices, LinearAlgebra, FFTW

julia> T = Toeplitz([0.;randn(499)], [0.;randn(999)]);

julia> x = randn(1000);

julia> length(T * x)
500

julia> Tfac = LinearAlgebra.factorize(T);

julia> Tfac * x
ERROR: MethodError: no method matching *(::ToeplitzMatrices.ToeplitzFactorization{Float64, Toeplitz{Float64, Vector{Float64}, Vector{Float64}}, ComplexF64, FFTW.cFFTWPlan{ComplexF64, -1, true, 1, UnitRange{Int64}}}, ::Vector{Float64})

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:578
  *(::StridedMatrix{T}, ::StridedVector{S}) where {T<:Union{Float32, Float64, ComplexF32, ComplexF64}, S<:Real}
   @ LinearAlgebra ~/julia-1.9.0-linux-x86_64/julia-1.9.0/share/julia/stdlib/v1.9/LinearAlgebra/src/matmul.jl:49
  *(::Union{SparseArrays.AbstractSparseMatrixCSC{TA, Ti}, SubArray{TA, 2, <:SparseArrays.AbstractSparseMatrixCSC{TA, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}, I}} where I<:AbstractUnitRange} where Ti, ::Union{StridedVector, BitVector}) where TA
   @ SparseArrays ~/julia-1.9.0-linux-x86_64/julia-1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/linalg.jl:50
  ...

Stacktrace:
 [1] top-level scope
   @ REPL[26]:1

I tried to look at how the library itself calls the mul! methods but it's weird (I don't understand, for example, what exactly alpha and beta are doing).

julia> methodswith(ToeplitzMatrices.ToeplitzFactorization)
[1] mul!(C::StridedMatrix{T} where T, A::ToeplitzMatrices.ToeplitzFactorization, B::StridedMatrix{T} where T, α::Number, β::Number) @ ToeplitzMatrices ~/.julia/packages/ToeplitzMatrices/5gW1c/src/linearalgebra.jl:88
[2] mul!(y::StridedVector{T} where T, A::ToeplitzMatrices.ToeplitzFactorization, x::StridedVector{T} where T, α::Number, β::Number) @ ToeplitzMatrices ~/.julia/packages/ToeplitzMatrices/5gW1c/src/linearalgebra.jl:42

It seems to work when passing 1.0 for alpha and beta, but it feels weird to use it without understand what's going on.

julia> result = zeros(Float64, 500); mul!(result, Tfac, x, 1., 1.,) == T*x
true

How does one use this? Maybe it could be documented in the README or elsewhere?

Another thing I'm confused about is whether the FFT functionality actually depends on something like FFTW being using'ed. To me it seems right now that adding ToeplitzMatrices to a project automatically includes FFTW through the DSP dependency. So while it may be possible to swap out FFTW for another library (I don't know if it is), it does not seem to be the case that a user needs to opt-in to the FFT functionality by making such a choice. Rather, FFTW is the default choice added as a dependency by ToeplitzMatrices and used from the beginning. If this is true, the README seems a bit unclear

To be able to use these methods, you have to install and load a package that implements the AbstractFFTs.jl  interface such as FFTW.jl

I read this as saying the FFT functionality is opt-in.

[adomasbaliuka]

I also found out that the function

mul!(y::StridedVector{T} where T, A::ToeplitzMatrices.ToeplitzFactorization, x::StridedVector{T} where T, α::Number, β::Number) @ ToeplitzMatrices ~/.julia/packages/ToeplitzMatrices/5gW1c/src/linearalgebra.jl:42

is not thread-safe (which of course nobody claimed it was).

[dlfivefifty]

I'm sure a PR to make it thread-safe would be merged.

[adomasbaliuka]

The mul! function is not thread-safe because it modifies the buffer inside the ToeplitzFactorization object. That means the matrix vector product T * x is thread-safe, since the factorization is computed per-multiplication. This also means that making it thread-safe has a runtime overhead and only makes a difference if using the mul!(::StridedVector{T} where T, ::ToeplitzFactorization, ::StridedVector{T} where T, ::Number, ::Number) method directly with a pre-computed factorization that is reused across threads.

The use-case is niche (and this mul! method is undocumented) so I'm rather worried to be wasting everyone's time by submitting a PR here. I'm also unsure what's the best way to make it thread-safe, as well as what other parts of the library are thread safe or not.

That said, here we go: #103