overlap now works well (adjusted maxevals)

src/_specfuncs.jl

@@ -65,35 +65,28 @@ end
 
 # Adapted from https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/polynomial/legendre.py#L832-L914
 function legval(x, c::AbstractVector)
-    # Pad the coefficient array if it contains only one (i.e. the constant 
+    # Pad the coefficient vector if it contains only one (i.e. the constant 
     # polynomial's) coefficient.
     length(c) ≥ 2 || (c = [c; zero(c)])
-
-    c0 = c[end - 1]
-    c1 = c[end]
-    tmp = copy(c0)
-
     nd = length(c)
-    @inbounds for i in 2:(length(c) - 1)
-        nd -= 1
-        invnd = inv(nd)
-        tmp = c0
-        c0 = c[end - i] - c1 * (1 - invnd)
-        c1 = tmp + c1 * x * (2 - invnd)
+
+    c0, c1 = c[nd - 1], c[nd]
+    @inbounds for j in (nd-2):-1:1
+        k = j / (j + 1)
+        c0, c1 = c[j] - c1 * k, c0 + c1 * x * (k + 1)
     end
-
     return c0 + c1 * x
 end
 
 # Adapted from https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/polynomial/legendre.py#L1126-L1176
 """
-    legvander(x, deg)
+legvander(x, deg)
 
 Pseudo-Vandermonde matrix of degree `deg`.
 """
 function legvander(x::Array{T,N}, deg::Integer) where {T,N}
     deg ≥ 0 || throw(DomainError(deg, "legvander needs a non-negative degree"))
-
+    
     # leading dimensions
     l = ntuple(_ -> :, N)
 

src/gauss.jl

@@ -1,4 +1,4 @@
-using QuadGK: gauss
+using QuadGK: gauss, kronrod
 using ._SpecFuncs: legvander
 
 export legendre, legendre_collocation, Rule, piecewise, quadrature, reseat
@@ -26,7 +26,8 @@ struct Rule{T}
         all(≤(b), x) || error("x must be ≤ b")
         all(≥(a), x) || error("x must be ≥ a")
         issorted(x) || error("x must be strictly increasing")
-        length(x) == length(w) || throw(DimensionMismatch("x and w must have the same length"))
+        length(x) == length(w) ||
+            throw(DimensionMismatch("x and w must have the same length"))
         return new{eltype(x)}(x, w, a, b)
     end
 end
@@ -112,3 +113,27 @@ end
 function Base.convert(::Type{Rule{T}}, rule::Rule{S}) where {T,S<:AbstractFloat}
     return Rule(T.(rule.x), T.(rule.w), T(rule.a), T(rule.b))
 end
+
+"""
+    NestedRule{T}
+
+Nested quadrature rule.
+"""
+struct NestedRule{T}
+    rule::Rule{T}
+    v::Vector{T}
+end
+
+function reseat(nrule::NestedRule, a, b)
+    newrule = reseat(nrule.rule, a, b)
+    newv = (b - a) / (rule.b - rule.a) * nrule.v
+    return NestedRule(newrule, newv)
+end
+
+function kronrod_21()
+    x, w, v = kronrod(10)
+    append!(x, -reverse(x)[2:end])
+    append!(w, reverse(w)[2:end])
+    append!(v, reverse(v))
+    return NestedRule(Rule(x, w), v)
+end

src/poly.jl

@@ -65,11 +65,11 @@ end
 (poly::PiecewiseLegendrePoly)(x) = _evaluate(poly, x)
 
 function _evaluate(poly::PiecewiseLegendrePoly, x::Number)
-    i, tildex= _split(poly, x)
-    return legval(tildex, view(poly.data,:, i)) * poly.norm[i]
+    i, tildex = _split(poly, x)
+    return legval(tildex, view(poly.data, :, i)) * poly.norm[i]
 end
 
-function _evaluate(poly::PiecewiseLegendrePoly, xs::Vector{T}) where T <: Real
+function _evaluate(poly::PiecewiseLegendrePoly, xs::AbstractVector)
     return poly.(xs)
 end
 
@@ -85,7 +85,7 @@ Given the function `f`, evaluate the integral::
 using adaptive Gauss-Legendre quadrature.
 """
 function overlap(poly::PiecewiseLegendrePoly, f; rtol=2.3e-16, return_error=false)
-    int_result, int_error = quadgk(x -> poly(x) * f(x), poly.knots...; rtol)
+    int_result, int_error = quadgk(x -> poly(x) * f(x), poly.knots...; rtol, order=10, maxevals=10^4)
     if return_error
         return int_result, int_error
     else
@@ -144,8 +144,8 @@ function roots(poly::PiecewiseLegendrePoly{T}; tol=1e-11) where {T}
     return rts
 end
 
-function check_domain(poly, x::Union{Number,Array})
-    all(poly.xmin ≤ x ≤ poly.xmax) || throw(DomainError(x, "x is outside the domain"))
+function check_domain(poly::PiecewiseLegendrePoly, x)
+    poly.xmin ≤ x ≤ poly.xmax || throw(DomainError(x, "x is outside the domain"))
     return true
 end
 

test/_specfuncs.jl

@@ -9,12 +9,10 @@ sphericalbesselj_M(n, x) = weval(W"SphericalBesselJ"(W"n", W"x"); n, x)
         @test !isnan(SparseIR.sphericalbesselj(4, 1e20))
     end
 
-    @testset "accuracy" begin
-        for pn in 0:7, px in 0:25
-            n = 2^pn
-            x = float(2^px)
-            @test SparseIR.sphericalbesselj(n, x) ≈ sphericalbesselj_M(n, x)
-            @test SparseIR.sphericalbesselj(n + 1, x) ≈ sphericalbesselj_M(n + 1, x)
-        end
+    @testset "accuracy (with pn = $pn, px = $px)" for pn in 0:7, px in 0:25
+        n = 2^pn
+        x = float(2^px)
+        @test SparseIR.sphericalbesselj(n, x) ≈ sphericalbesselj_M(n, x)
+        @test SparseIR.sphericalbesselj(n + 1, x) ≈ sphericalbesselj_M(n + 1, x)
     end
 end

test/augment.jl

@@ -3,55 +3,51 @@ using SparseIR
 using AssociatedLegendrePolynomials: Plm
 
 @testset "augment.jl" begin
-    @testset "MatsubaraConstBasis" begin
+    @testset "MatsubaraConstBasis with stat = $stat" for stat in (fermion, boson)
         beta = 2.0
         value = 1.1
-        for statistics in (fermion, boson)
-            b = MatsubaraConstBasis(statistics, beta; value=value)
-            shift::Int = Dict(fermion => 1, boson => 0)[statistics]
-            n = 2 .* collect(1:10) .+ shift
-            @test b.uhat(n) ≈ fill(value, 1, length(n))
-        end
+        b = MatsubaraConstBasis(stat, beta; value=value)
+        shift::Int = Dict(fermion => 1, boson => 0)[stat]
+        n = 2 .* collect(1:10) .+ shift
+        @test b.uhat(n) ≈ fill(value, 1, length(n))
     end
 
-    @testset "LegendreBasis" begin
-        for stat in (fermion, boson)
-            β = 1.0
-            Nl = 10
-            cl = sqrt.(2 * collect(0:(Nl - 1)) .+ 1)
-            basis = SparseIR.LegendreBasis(stat, β, Nl; cl=cl)
-            @test size(basis) == (Nl,)
-
-            τ = Float64[0, 0.1 * β, 0.4 * β, β]
-            uval = basis.u(τ)
-
-            ref = Matrix{Float64}(undef, Nl, length(τ))
-            for l in 0:(Nl - 1)
-                x = 2τ / β .- 1
-                ref[l + 1, :] .= cl[l + 1] * (√(2l + 1) / β) * Plm.(l, 0, x)
-            end
-            @test uval ≈ ref
-
-            sign = stat == fermion ? -1 : 1
-
-            # G(iv) = 1/(iv-pole)
-            # G(τ) = -e^{-τ*pole}/(1 + e^{-β*pole}) [F]
-            #        = -e^{-τ*pole}/(1 - e^{-β*pole}) [B]
-            pole = 1.0
-            τ_smpl = TauSampling(basis)
-            gτ = -exp.(-τ_smpl.sampling_points * pole) / (1 - sign * exp(-β * pole))
-            gl_from_τ = fit(τ_smpl, gτ)
-
-            matsu_smpl = MatsubaraSampling(basis)
-            giv = 1 ./ ((im * π / β) * matsu_smpl.sampling_points .- pole)
-            gl_from_matsu = fit(matsu_smpl, giv)
-
-            #println(maximum(abs.(gl_from_τ-gl_from_matsu)))
-            #println(maximum(abs.(gl_from_τ)))
-            #println("gl_from_τ", gl_from_τ[1:4])
-            #println("gl_from_matsu", gl_from_matsu[1:4])
-            @test isapprox(gl_from_τ, gl_from_matsu;
-                           atol=1e-10 * maximum(abs, gl_from_matsu), rtol=0)
+    @testset "LegendreBasis with stat = $stat" for stat in (fermion, boson)
+        β = 1.0
+        Nl = 10
+        cl = sqrt.(2 * collect(0:(Nl - 1)) .+ 1)
+        basis = SparseIR.LegendreBasis(stat, β, Nl; cl=cl)
+        @test size(basis) == (Nl,)
+
+        τ = Float64[0, 0.1 * β, 0.4 * β, β]
+        uval = basis.u(τ)
+
+        ref = Matrix{Float64}(undef, Nl, length(τ))
+        for l in 0:(Nl - 1)
+            x = 2τ / β .- 1
+            ref[l + 1, :] .= cl[l + 1] * (√(2l + 1) / β) * Plm.(l, 0, x)
         end
+        @test uval ≈ ref
+
+        sign = stat == fermion ? -1 : 1
+
+        # G(iv) = 1/(iv-pole)
+        # G(τ) = -e^{-τ*pole}/(1 + e^{-β*pole}) [F]
+        #        = -e^{-τ*pole}/(1 - e^{-β*pole}) [B]
+        pole = 1.0
+        τ_smpl = TauSampling(basis)
+        gτ = -exp.(-τ_smpl.sampling_points * pole) / (1 - sign * exp(-β * pole))
+        gl_from_τ = fit(τ_smpl, gτ)
+
+        matsu_smpl = MatsubaraSampling(basis)
+        giv = 1 ./ ((im * π / β) * matsu_smpl.sampling_points .- pole)
+        gl_from_matsu = fit(matsu_smpl, giv)
+
+        #println(maximum(abs.(gl_from_τ-gl_from_matsu)))
+        #println(maximum(abs.(gl_from_τ)))
+        #println("gl_from_τ", gl_from_τ[1:4])
+        #println("gl_from_matsu", gl_from_matsu[1:4])
+        @test isapprox(gl_from_τ, gl_from_matsu;
+                       atol=1e-10 * maximum(abs, gl_from_matsu), rtol=0)
     end
-end
+end

test/kernel.jl

@@ -2,32 +2,31 @@ using Test
 using SparseIR
 
 @testset "kernel.jl" begin
-    @testset "accuracy" begin
-        for K in [LogisticKernel(9.0),
-                  RegularizedBoseKernel(8.0),
-                  LogisticKernel(120_000.0),
-                  RegularizedBoseKernel(127_500.0),
-                  get_symmetrized(LogisticKernel(40_000.0), -1),
-                  get_symmetrized(RegularizedBoseKernel(35_000.0), -1)]
-            T = Float32
-            T_x = Float64
+    @testset "accuracy with K = $K" for K in (LogisticKernel(9.0),
+                                              RegularizedBoseKernel(8.0),
+                                              LogisticKernel(120_000.0),
+                                              RegularizedBoseKernel(127_500.0),
+                                              get_symmetrized(LogisticKernel(40_000.0), -1),
+                                              get_symmetrized(RegularizedBoseKernel(35_000.0),
+                                                              -1))
+        T = Float32
+        T_x = Float64
 
-            rule = convert(Rule{T}, legendre(10))
-            hints = sve_hints(K, 2.2e-16)
-            gauss_x = piecewise(rule, segments_x(hints))
-            gauss_y = piecewise(rule, segments_y(hints))
-            ϵ = eps(T)
-            tiny = floatmin(T) / ϵ
+        rule = convert(Rule{T}, legendre(10))
+        hints = sve_hints(K, 2.2e-16)
+        gauss_x = piecewise(rule, segments_x(hints))
+        gauss_y = piecewise(rule, segments_y(hints))
+        ϵ = eps(T)
+        tiny = floatmin(T) / ϵ
 
-            result = matrix_from_gauss(K, gauss_x, gauss_y)
-            result_x = matrix_from_gauss(K,
-                                         convert(Rule{T_x}, gauss_x),
-                                         convert(Rule{T_x}, gauss_y))
-            magn = maximum(abs, result_x)
+        result = matrix_from_gauss(K, gauss_x, gauss_y)
+        result_x = matrix_from_gauss(K,
+                                     convert(Rule{T_x}, gauss_x),
+                                     convert(Rule{T_x}, gauss_y))
+        magn = maximum(abs, result_x)
 
-            @test result ≈ result_x atol = 2magn * ϵ rtol = 0
-            reldiff = @. ifelse(abs(result) < tiny, 1, result / result_x)
-            @test all(isapprox.(reldiff, 1, atol=100ϵ, rtol=0))
-        end
+        @test result ≈ result_x atol = 2magn * ϵ rtol = 0
+        reldiff = @. ifelse(abs(result) < tiny, 1, result / result_x)
+        @test all(isapprox.(reldiff, 1, atol=100ϵ, rtol=0))
     end
 end

test/matsubara.jl

@@ -1,39 +1,38 @@
 using Test
 using SparseIR
 
-const ε = (@isdefined Double64) ? nothing : 1e-15
-
 @testset "matsubara.jl" begin
-    @testset "single pole" begin
-        for stat in (fermion, boson), Λ in (1e1, 1e4)
-            wmax = 1.0
-            pole = 0.1 * wmax
-            β = Λ / wmax
-            basis = FiniteTempBasis(stat, β, wmax, ε; sve_result=sve_logistic[Λ])
+    @testset "single pole with stat = $stat, Λ = $Λ" for stat in (fermion, boson),
+                                                         Λ in (1e1, 1e4)
+
+        ε = (@isdefined Double64) ? nothing : 1e-15
+        wmax = 1.0
+        pole = 0.1 * wmax
+        β = Λ / wmax
+        basis = FiniteTempBasis(stat, β, wmax, ε; sve_result=sve_logistic[Λ])
 
-            stat_shift = (stat == fermion) ? 1 : 0
-            weight = (stat == fermion) ? 1 : 1 / tanh(0.5 * Λ * pole / wmax)
-            gl = -basis.s .* basis.v(pole) * weight
+        stat_shift = (stat == fermion) ? 1 : 0
+        weight = (stat == fermion) ? 1 : 1 / tanh(0.5 * Λ * pole / wmax)
+        gl = -basis.s .* basis.v(pole) * weight
 
-            func_G(n) = 1 / (im * (2n + stat_shift) * π / β - pole)
+        func_G(n) = 1 / (im * (2n + stat_shift) * π / β - pole)
 
-            # Compute G(iwn) using unl
-            matsu_test = Int[-1, 0, 1, 1e2, 1e4, 1e6, 1e8, 1e10, 1e12]
-            prj_w = transpose(basis.uhat(2matsu_test .+ stat_shift))
-            Giwn_t = prj_w * gl
+        # Compute G(iwn) using unl
+        matsu_test = Int[-1, 0, 1, 1e2, 1e4, 1e6, 1e8, 1e10, 1e12]
+        prj_w = transpose(basis.uhat(2matsu_test .+ stat_shift))
+        Giwn_t = prj_w * gl
 
-            # Compute G(iwn) from analytic expression
-            Giwn_ref = func_G.(matsu_test)
+        # Compute G(iwn) from analytic expression
+        Giwn_ref = func_G.(matsu_test)
 
-            magnitude = maximum(abs, Giwn_ref)
-            diff = abs.(Giwn_t - Giwn_ref)
-            tol = max(10 * last(basis.s) / first(basis.s), 1e-10)
+        magnitude = maximum(abs, Giwn_ref)
+        diff = abs.(Giwn_t - Giwn_ref)
+        tol = max(10 * last(basis.s) / first(basis.s), 1e-10)
 
-            # Absolute error
-            @test maximum(diff ./ magnitude) < tol
+        # Absolute error
+        @test maximum(diff ./ magnitude) < tol
 
-            # Relative error
-            @test maximum(abs, diff ./ Giwn_ref) < tol
-        end
+        # Relative error
+        @test maximum(abs, diff ./ Giwn_ref) < tol
     end
 end