From d6bb4f2f549d6bde06b7a8c14999491b25337cc6 Mon Sep 17 00:00:00 2001 From: Haroun Meghaichi Date: Thu, 26 Jun 2025 18:46:05 -0400 Subject: [PATCH 1/6] re-use PP1 and PP2 to reduce allocations --- src/gausslegendre.jl | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/gausslegendre.jl b/src/gausslegendre.jl index 4049b51..ba22061 100644 --- a/src/gausslegendre.jl +++ b/src/gausslegendre.jl @@ -207,14 +207,12 @@ function rec(n) x0 = asy(n)[1] x = x0[1:n ÷ 2 + 1] # 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 +230,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 +249,10 @@ 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 \ No newline at end of file From 07d26b4363064586d88ba2f87ee5faf42d2b654c Mon Sep 17 00:00:00 2001 From: Haroun Meghaichi Date: Fri, 27 Jun 2025 16:12:51 -0400 Subject: [PATCH 2/6] add leg_initial_guess: Basically the same as asy(n) but with less allocations --- src/gausslegendre.jl | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/src/gausslegendre.jl b/src/gausslegendre.jl index ba22061..b12cb95 100644 --- a/src/gausslegendre.jl +++ b/src/gausslegendre.jl @@ -203,9 +203,10 @@ 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=similar(x),similar(x) @@ -255,4 +256,16 @@ end @inbounds @simd for i in eachindex(x) x[i] -= PP1[i] / PP2[i] end +end + +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 + 1) >> 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 From dc3e83a07c484db1da7e184a081b835c3df2037c Mon Sep 17 00:00:00 2001 From: Haroun Meghaichi Date: Fri, 27 Jun 2025 16:22:06 -0400 Subject: [PATCH 3/6] the size of the initial guess is n\div 2 +1 --- src/gausslegendre.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gausslegendre.jl b/src/gausslegendre.jl index b12cb95..5a56749 100644 --- a/src/gausslegendre.jl +++ b/src/gausslegendre.jl @@ -262,7 +262,7 @@ 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 + 1) >> 1 + m = (n ÷2) +1 a = besselZeroRoots(m) rmul!(a, 1 / (n + 0.5)) x = legpts_nodes(n, a) From 0bd71f1eeae1af9fdf6db5b893f24b177b0b4b36 Mon Sep 17 00:00:00 2001 From: Haroun Meghaichi Date: Fri, 27 Jun 2025 16:38:15 -0400 Subject: [PATCH 4/6] add inline flag to initial guess function. --- src/gausslegendre.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gausslegendre.jl b/src/gausslegendre.jl index 5a56749..d8d9457 100644 --- a/src/gausslegendre.jl +++ b/src/gausslegendre.jl @@ -258,7 +258,7 @@ end end end -function leg_initial_guess(n) +@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. From 52e07c163a73c94dd30424cdf72201525544afff Mon Sep 17 00:00:00 2001 From: Haroun Meghaichi Date: Fri, 27 Jun 2025 17:00:28 -0400 Subject: [PATCH 5/6] CI.yml: Upgrade actions/cachev1 (deprecated) to v3 --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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: From 960629942b20dc671974b80f45d9db9346183830 Mon Sep 17 00:00:00 2001 From: Haroun Meghaichi Date: Fri, 27 Jun 2025 17:42:30 -0400 Subject: [PATCH 6/6] Aqua keyword arg: check_extras=false. This is a temporary fix to the errors in the CI since Random and Test have no compat entries --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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")