diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 55bffd9..10c87e5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -30,7 +30,7 @@ jobs: with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: actions/cache@v1 + - uses: actions/cache@v3 env: cache-name: cache-artifacts with: diff --git a/src/gausslegendre.jl b/src/gausslegendre.jl index 4049b51..d8d9457 100644 --- a/src/gausslegendre.jl +++ b/src/gausslegendre.jl @@ -203,18 +203,17 @@ end function rec(n) # COMPUTE GAUSS-LEGENDRE NODES AND WEIGHTS USING NEWTON'S METHOD. # THREE-TERM RECURENCE IS USED FOR EVALUATION. COMPLEXITY O(n^2). - # Initial guesses: - x0 = asy(n)[1] - x = x0[1:n ÷ 2 + 1] + + # Initial guess: + x=leg_initial_guess(n) + # Perform Newton to find zeros of Legendre polynomial: - PP1, PP2 = innerRec(n, x) - @inbounds @simd for i in 1:length(x) - x[i] -= PP1[i] / PP2[i] - end - # One more Newton for derivatives: - PP1, PP2 = innerRec(n, x) - @inbounds @simd for i in 1:length(x) - x[i] -= PP1[i] / PP2[i] + PP1,PP2=similar(x),similar(x) + + # Two iterations of Newton's method + for _=1:2 + innerRec!(PP1,PP2,n, x) + newt_step!(x,PP1,PP2) end # Use symmetry to get the other Legendre nodes and weights: @@ -232,11 +231,9 @@ function rec(n) return x, PP2 end -function innerRec(n, x) +@inline function innerRec!(myPm1,myPPm1,n, x) # EVALUATE LEGENDRE AND ITS DERIVATIVE USING THREE-TERM RECURRENCE RELATION. N = size(x, 1) - myPm1 = Array{Float64}(undef, N) - myPPm1 = Array{Float64}(undef, N) @inbounds for j = 1:N xj = x[j] Pm2 = 1.0 @@ -253,3 +250,22 @@ function innerRec(n, x) end return myPm1, myPPm1 end + +@inline function newt_step!(x,PP1,PP2) + # In-place iteration of Newton's method. + @inbounds @simd for i in eachindex(x) + x[i] -= PP1[i] / PP2[i] + end +end + +@inline function leg_initial_guess(n) + # Returns an approximation of the first n÷2+1 roots of the Legendre polynomial. + # The following is equivalent to "x0=asy(n);x = x0[1:n ÷ 2 + 1]" but it avoids unnecessary calculations. + + m = (n ÷2) +1 + a = besselZeroRoots(m) + rmul!(a, 1 / (n + 0.5)) + x = legpts_nodes(n, a) + rmul!(x,-1.0) + return x +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 0b60392..6409fc5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,7 @@ using FastGaussQuadrature using Test, Aqua, LinearAlgebra, Random, SpecialFunctions -Aqua.test_all(FastGaussQuadrature) +Aqua.test_all(FastGaussQuadrature;deps_compat = (; check_extras=false)) @testset "FastGaussQuadrature.jl" begin include("test_gausschebyshev.jl")