Linear System Solver for Toeplitz slow and does not agree with LU-algorithm

[Firionus]

I used your package for the first time today and experienced strange results solving linear systems of Toeplitz matrices:

julia> using ToeplitzMatrices, BenchmarkTools, LinearAlgebra

julia> x = rand(1_000); y=rand(1_000); z=rand(1_000);

julia> x[1]=y[1]
0.8516822317176573

julia> A = Toeplitz(x,y);

julia> @benchmark A\z
BenchmarkTools.Trial:
  memory estimate:  306.00 MiB
  allocs estimate:  2193266
  --------------
  minimum time:     6.468 s (1.05% GC)
  median time:      6.468 s (1.05% GC)
  mean time:        6.468 s (1.05% GC)
  maximum time:     6.468 s (1.05% GC)
  --------------
  samples:          1
  evals/sample:     1

julia> lu_factored = lu(A);

julia> issuccess(lu_factored)
true

julia> @benchmark lu(A)\z
BenchmarkTools.Trial:
  memory estimate:  7.65 MiB
  allocs estimate:  5
  --------------
  minimum time:     31.134 ms (0.00% GC)
  median time:      36.463 ms (0.00% GC)
  mean time:        37.242 ms (3.35% GC)
  maximum time:     97.156 ms (56.64% GC)
  --------------
  samples:          135
  evals/sample:     1

Not only does the standard way of solving the Toeplitz matrix take very long, it is also different from the LU-factored solution, as you can see in this plot of the solutions:

Plot was created by:

julia> direct_solution = A\z;

julia> using PyPlot;

julia> plot(direct_solution); plot(lu_factored\z); legend(["direct", "LU"]);

(v1.2) pkg> status ToeplitzMatrices
    Status `C:\Users\Johannes\.julia\environments\v1.2\Project.toml`
  [7a1cc6ca] FFTW v1.1.0
  [c751599d] ToeplitzMatrices v0.6.1

[andreasnoack]

The default solver for Toeplitz matrices is an iterative solver since it's generally not possible to benefit from the Toeplitz structure in a direct solver (to my knowledge). The iterative solver is pretty slow for smaller problems, but you can try with n=16_000. The iterative solver is also not as accurate as a direct solver so if you can afford it, you should probably use LU.