Tests fail with more accurate sinc differentiation matrices

[MikaelSlevinsky]

Prior to today, the second order differentiation matrix's infinity-norms of the error on the sinc points was:

julia> using SincFun # checkout commit 4a656dc29a34247c65388d11d513a5982b2f256d

julia> k=collect(-500.0:500.0);

julia> norm(sinc(2,k)-convert(Vector{Float64},sinc(2,convert(Vector{BigFloat},k))),Inf)
1.734723475976807e-18

julia> norm(sinc(2,k)./convert(Vector{Float64},sinc(2,convert(Vector{BigFloat},k)))-1,Inf)
3.3306690738754696e-16

Now, I updated the code to make it essentially allocation-free, and the infinity-norm is also smaller:

julia> using SincFun # checkout master

julia> k=collect(-500.0:500.0);

julia> norm(sinc(2,k)-convert(Vector{Float64},sinc(2,convert(Vector{BigFloat},k))),Inf)
4.336808689942018e-19

julia> norm(sinc(2,k)./convert(Vector{Float64},sinc(2,convert(Vector{BigFloat},k)))-1,Inf)
2.220446049250313e-16

The tests fail because it's not clear how many eigenvalues will be returned, just the ones below the threshold tolerance. And somehow the results are worse in that there are fewer eigenvalues returned than before?

[pjgaudre]

The DESCM on semi-infinite domains where the eigenfunctions decays algebraically at 0 and exponentially at infinity have badly conditioned matrices because the condition number of the matrix D2 exceeds 10^16 quite rapidly. Since LAPACK for Julia only deals with double precision, it makes sense that we would not have that many eigenvalues that converged. However, it is strange that the method is doing worse than before.

[MikaelSlevinsky]

Can we use our own Arnoldi or Lanczos? It should be straightforward to use SincFun's SincMatrix to get fast and data-sparse matrix-vector multiply.

[pjgaudre]

I tried coding my own version of a generalized eigenvalue solver about a year ago and found it incredibly challenging. Eigenvalue solvers are very finicky it's no easy task. But if you're up to the task, that would certainly be useful!

[MikaelSlevinsky]

SincFun has Lanczos and steig! (symmetric tridiagonal) methods. I could help by increasing their performance, but you can only get so much out of eigvals.

[pjgaudre]

I agree. If you can get your Lanczos algorithm to be just as performant or better, we can definitely make the switch.