inplace versions

Author: dlfivefifty

src/MultivariateSingularIntegrals.jl

@@ -68,11 +68,10 @@ end
 
 
 
-function rec_rhs_1(z, m, n)
+function rec_rhs_1!(F::AbstractMatrix{T}, z) where T
+    m,n = size(F)
     x,y = reim(z)
-    T= complex(float(typeof(z)))
     πT = convert(T, π)
-    F = zeros(T,m,n)
     if x < -1 && -1 < y < 1
         C_y = Ultraspherical{T}(-1/2)[y,2:n+1]
         F[1,:] .-=  (4im*πT * x) .* C_y
@@ -97,9 +96,9 @@ end
 
 
 
-function rec_rhs_2(z, m, n)
+function rec_rhs_2!(F::AbstractMatrix{T}, z) where T
+    m,n = size(F)
     x,y = reim(z)
-    T= T= complex(float(typeof(z)))
     πT = convert(T, π)
     if -1 < x < 1 && -1 ≤ y < 1
         C_x = Ultraspherical{T}(-3/2)[x,3:m+2]
@@ -154,17 +153,26 @@ The bottom right of the returned matrix is zero. For a square truncation compute
 """
 
 function logkernelsquare(z, n)
-    F_1 = rec_rhs_1(z, n, n)
-    F_2 = rec_rhs_2(z, n, n)
-    A = zeros(complex(float(eltype(z))),n,n)
+    T = complex(float(eltype(z)))
+    logkernelsquare!(zeros(T,n,n), z, zeros(T,n,n), zeros(T,n,n))
+end
+
+
+function logkernelsquare!(A::AbstractMatrix{T}, z, F_1, F_2) where T
+    m,n = size(A)
+    @assert m  == n
+    rec_rhs_1!(F_1, z)
+    rec_rhs_2!(F_2, z)
     A[1,1] = L₀₀(z)
     logkernelsquare_populatefirstcolumn!(A, z, F_1, F_2)
     logkernelsquare_populatefirstrow!(A, z, F_1, F_2)
-    F = F_2 - F_1
+
+    F = F_1 # reuse the memory
+    F .= F_2 .- F_1
 
     # 2nd row/column
 
-    for k = 1:n-2
+    for k = 1:m-2
         A[k+1,2] = im*(F[k+1,1] + (A[k,1] - A[k+2,1])/(2k+1))
     end
 

test/runtests.jl

@@ -1,13 +1,30 @@
 using MultivariateSingularIntegrals, QuadGK, ClassicalOrthogonalPolynomials, Test
 
+@testset "accurate for low order" begin
+    Z = (2, 2im, 2+2im, 2-2im, -2 - 2im, -1.001-1.5im, -0.999-1.5im, -0.5-2im, -1-1.5im, -2+2im, -2, 0.1+0.2im, 0.1+im, 0.1-im, 1+0.1im, -1+0.1im, 1+im, -1+im, -1-im, 1-im)
+    n = 6
+    for z in Z
+        L = logkernelsquare(z,n)
+        for j = 0:n-1, k=0:n-j-1
+            q = quadgk(s -> quadgk(t -> iszero(s+im*t) ? 0.0+0.0im : log(z-(s+im*t)) * legendrep(k,s) * legendrep(j,t), -1, 1, atol=1E-12)[1], -1, 1, atol=1E-12)[1]
+            @test q ≈ L[k+1,j+1] atol=1E-10
+        end
+    end
+end
 
-Z = (2, 2im, 2+2im, 2-2im, -2 - 2im, -1.001-1.5im, -0.999-1.5im, -0.5-2im, -1-1.5im, -2+2im, -2, 0.1+0.2im, 0.1+im, 0.1-im, 1+0.1im, -1+0.1im, 1+im, -1+im, -1-im, 1-im)
-n = 6
-for z in Z
-    L = logkernelsquare(z,n)
-    for j = 0:n-1, k=0:n-j-1
-        q = quadgk(s -> quadgk(t -> iszero(s+im*t) ? 0.0+0.0im : log(z-(s+im*t)) * legendrep(k,s) * legendrep(j,t), -1, 1, atol=1E-12)[1], -1, 1, atol=1E-12)[1]
-        @test q ≈ L[k+1,j+1] atol=1E-10
+@testset "BigFloat" begin
+    setprecision(2000) do
+        for z in (big(2.0), big(2.0im), big(2.0+2im), big(2.0-2im), big(-2.0 - 2im))
+            P = Legendre{BigFloat}()
+            n = 100
+            M = Diagonal((P'P)[1:n,1:n])
+            @test P[big(0.), 1:n]' *  inv(M)  * logkernelsquare(z, n) * inv(M) * P[big(0.), 1:n] ≈ log(z) atol=1E-40
+        end
     end
 end
 
+n = 4000
+z = 2.0
+@profview logkernelsquare(z, n)
+
+x = range(-1,1,n); @time log.(abs.(x .+ im .* x'));