PEPs in Chebyshev basis

[jarlebring]

This adds representation of a PEP in a Chebyshev basis (scaled to arbitrary interval). It supports creation via interpolation of another NEP in Chebyshev nodes, and an implementation of polyeig which uses a Chebyshev companion linearization:

julia> nep=nep_gallery("dep0",5);
julia> n=size(nep,1);
julia> chebpep=ChebPEP(nep,13,-1,1); # Create ChebPEP based on interpolation in chebyshev nodes
julia> compute_Mder(nep,0)
5×5 Array{Float64,2}:
  0.790529   1.72038   -0.182033   0.293983  -0.302562
  0.47053    0.473334  -0.819734  -1.9292     0.66847 
  0.120707   0.479604  -1.37354   -2.23445    0.979773
  0.16538   -0.810488   2.40225    2.09744   -0.155803
 -0.17606    2.26544    1.05071   -1.25819    0.856748
julia> compute_Mder(chebpep,0)
5×5 Array{Float64,2}:
  0.790529   1.72038   -0.182033   0.293983  -0.302562
  0.47053    0.473334  -0.819734  -1.9292     0.66847 
  0.120707   0.479604  -1.37354   -2.23445    0.979773
  0.16538   -0.810488   2.40225    2.09744   -0.155803
 -0.17606    2.26544    1.05071   -1.25819    0.856748
julia> v0=ones(n);
julia> (λ2,V2)=iar(chebpep,neigs=2,errmeasure=ResidualErrmeasure,v=v0)
(Complex{Float64}[-0.358719+4.25641e-17im, 0.834735-2.3048e-16im], Complex{Float64}[-0.297881+2.63211e-16im 0.534069+2.23465e-16im; 0.149929+3.37404e-17im 0.0702221+1.40189e-16im; … ; 0.524012+1.09274e-17im 0.238219-7.1747e-17im; 0.641947+6.04079e-17im 0.420887-1.79986e-16im], Complex{Float64}[0.447214+0.0im 0.721426-0.0im … -0.0108986+1.07354e-11im 0.0258408-9.68579e-11im; 0.447214+0.0im -0.250357-0.0im … 0.0279483-2.95601e-11im -0.00884929+5.71576e-12im; … ; 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im])
julia> norm(compute_Mlincomb(nep,λ2[1],V2[:,1]))
5.403860318376269e-14
julia> (λ,V)=polyeig(chebpep);  # polyeig with chebyshev basis
julia> j=argmin(abs.(λ .- λ2[1])); # See if we found what we found with iar
julia> λ[j]
-0.35871894596860465 + 0.0im
julia> λ2[1]
-0.35871894596860465 + 4.256411657408781e-17im

Still remaining:

• polyeig(ChebPEP) should not only support the unit interval.
• There are two ways to compute chebyshev polys: the monomial formula or the cosine-formula. If we only consider efficiency, the best would probably be to make a cut-off and switch to cosine-formula for the higher poly's.
• *fixed by new acos(LowerTriangular) * Derivatives currently do not work since acos(A) rounds in a strange way (cf julia issue 32721)

julia> compute_Mder(chebpep,0.1,1);
julia> compute_Mder(chebpep,0.1,2);
ERROR: InexactError: Float64(0.0 - 5.551115123125783e-16im)