Mercurial > repos > public > sbplib_julia

changeset 1554:05e0bcf7d7ca refactor/equidistant_grid/signature

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 13 Apr 2024 23:46:06 +0200
parents 43aaf710463e (current diff) f1eacb923f45 (diff)
children 7e165cc0eb68
files
diffstat 5 files changed, 128 insertions(+), 35 deletions(-) [+]
line wrap: on
line diff
--- a/Notes.md	Thu Apr 11 22:31:04 2024 +0200
+++ b/Notes.md	Sat Apr 13 23:46:06 2024 +0200
@@ -348,3 +348,11 @@
 * This would impact how tensor grid works.
 * To make things homogenous maybe these index collections should be used for the more simple grids too.
 * The function `to_indices` from Base could be useful to implement for `IndexCollection`
+
+
+## Stencil application pipeline
+We should make sure that `@inbounds` and `Base.@propagate_inbounds` are
+applied correctly throughout the stack. When testing the performance of
+stencil application on the bugfix/sbp_operators/stencil_return_type branch
+there seemed to be some strange results where such errors could be the
+culprit.
--- a/benchmark/Manifest.toml	Thu Apr 11 22:31:04 2024 +0200
+++ b/benchmark/Manifest.toml	Sat Apr 13 23:46:06 2024 +0200
@@ -1,19 +1,19 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.10.2"
 manifest_format = "2.0"
 project_hash = "25bba7b4a00465d5a2b00b589eb10e3301c31f2a"
 
 [[deps.AbstractTrees]]
-git-tree-sha1 = "faa260e4cb5aba097a73fab382dd4b5819d8ec8c"
+git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
 uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
-version = "0.4.4"
+version = "0.4.5"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "0fb305e0253fd4e833d486914367a2ee2c2e78d0"
+git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.0.1"
+version = "4.0.4"
 weakdeps = ["StaticArrays"]
 
     [deps.Adapt.extensions]
@@ -24,16 +24,18 @@
 version = "1.1.1"
 
 [[deps.ArrayInterface]]
-deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "bbec08a37f8722786d87bedf84eae19c020c4efa"
+deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "44691067188f6bd1b2289552a23e4b7572f4528d"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.7.0"
+version = "7.9.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
     ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
     ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceChainRulesExt = "ChainRules"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    ArrayInterfaceReverseDiffExt = "ReverseDiff"
     ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
     ArrayInterfaceTrackerExt = "Tracker"
 
@@ -41,7 +43,9 @@
     BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
     BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 
@@ -53,15 +57,15 @@
 
 [[deps.BenchmarkTools]]
 deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
-git-tree-sha1 = "f1f03a9fa24271160ed7e73051fba3c1a759b53f"
+git-tree-sha1 = "f1dff6729bc61f4d49e140da1af55dcd1ac97b2f"
 uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.4.0"
+version = "1.5.0"
 
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
-git-tree-sha1 = "75bd5b6fc5089df449b5d35fa501c846c9b6549b"
+git-tree-sha1 = "c955881e3c981181362ae4088b35995446298b80"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.12.0"
+version = "4.14.0"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -70,7 +74,7 @@
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.1.0+0"
 
 [[deps.DataAPI]]
 git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe"
@@ -192,7 +196,7 @@
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.23+4"
 
 [[deps.OrderedCollections]]
 git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
@@ -218,15 +222,15 @@
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e"
+git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.4.1"
+version = "1.4.3"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -279,9 +283,9 @@
 
 [[deps.Static]]
 deps = ["IfElse"]
-git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc"
+git-tree-sha1 = "d2fdac9ff3906e27f7a618d47b676941baa6c80c"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.8.8"
+version = "0.8.10"
 
 [[deps.StaticArrayInterface]]
 deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"]
@@ -296,9 +300,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "f68dd04d131d9a8a8eb836173ee8f105c360b0c5"
+git-tree-sha1 = "bf074c045d3d5ffd956fa0a461da38a44685d6b2"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.1"
+version = "1.9.3"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/benchmark/benchmark_utils.jl	Thu Apr 11 22:31:04 2024 +0200
+++ b/benchmark/benchmark_utils.jl	Sat Apr 13 23:46:06 2024 +0200
@@ -23,7 +23,7 @@
 For control over what happens to the benchmark result datastructure see the
 different methods of [`run_benchmark`](@ref)
 """
-function main(;rev=nothing, target=nothing, baseline=nothing , kwargs...)
+function main(;rev=nothing, target=nothing, baseline=nothing, name=nothing, kwargs...)
     if !isnothing(rev)
         r = run_benchmark(rev; kwargs...)
     elseif !isnothing(baseline)
@@ -37,7 +37,7 @@
         r = run_benchmark(;kwargs...)
     end
 
-    file_path = write_result_html(r)
+    file_path = write_result_html(r; name)
     open_in_default_browser(file_path)
 end
 
@@ -137,9 +137,14 @@
     Mustache.render(io, template, Dict("title"=>dt, "content"=>content))
 end
 
-function write_result_html(r)
+function write_result_html(r; name=nothing)
     dt = Dates.format(PkgBenchmark.date(r), "yyyy-mm-dd HHMMSS")
-    file_path = joinpath(results_dir, dt*".html")
+
+    if isnothing(name)
+        file_path = joinpath(results_dir, dt*".html")
+    else
+        file_path = joinpath(results_dir, dt*" "*name*".html")
+    end
 
     open(file_path, "w") do io
         write_result_html(io, r)
--- a/src/SbpOperators/stencil.jl	Thu Apr 11 22:31:04 2024 +0200
+++ b/src/SbpOperators/stencil.jl	Sat Apr 13 23:46:06 2024 +0200
@@ -69,22 +69,29 @@
 end
 
 Base.@propagate_inbounds @inline function apply_stencil(s::Stencil, v::AbstractVector, i::Int)
-    w = zero(promote_type(eltype(s),eltype(v)))
-    @simd for k ∈ 1:length(s)
-        w += s.weights[k]*v[i + s.range[k]]
+    return sum(enumerate(s.weights)) do (k,w) #TBD: Which optimizations are needed here?
+        w*v[i + @inbounds s.range[k]]
     end
-
-    return w
 end
 
 Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil, v::AbstractVector, i::Int)
-    w = zero(promote_type(eltype(s),eltype(v)))
-    @simd for k ∈ length(s):-1:1
-        w += s.weights[k]*v[i - s.range[k]]
+    return sum(enumerate(s.weights)) do (k,w) #TBD: Which optimizations are needed here?
+        w*v[i - @inbounds s.range[k]]
     end
-    return w
 end
 
+# There are many options for the implementation of `apply_stencil` and
+# `apply_stencil_backwards`. Some alternatives were tried on the branch
+# bugfix/sbp_operators/stencil_return_type and can be found at the following
+# revision:
+#
+# * 237b980ffb91 (baseline)
+# * a72bab15228e (mapreduce)
+# * ffd735354d54 (multiplication)
+# * b5abd5191f2c (promote_op)
+# * 8d56846185fc (return_type)
+#
+
 function left_pad(s::Stencil, N)
     weights = LazyTensors.left_pad_tuple(s.weights, zero(eltype(s)), N)
     range = (first(s.range) - (N - length(s.weights))):last(s.range)
--- a/test/SbpOperators/stencil_test.jl	Thu Apr 11 22:31:04 2024 +0200
+++ b/test/SbpOperators/stencil_test.jl	Sat Apr 13 23:46:06 2024 +0200
@@ -1,8 +1,10 @@
 using Test
 using Sbplib.SbpOperators
+using StaticArrays
 import Sbplib.SbpOperators.Stencil
 import Sbplib.SbpOperators.NestedStencil
 import Sbplib.SbpOperators.scale
+import Sbplib.SbpOperators: apply_stencil, apply_stencil_backwards
 
 @testset "Stencil" begin
     s = Stencil(-2:2, (1.,2.,2.,3.,4.))
@@ -44,6 +46,52 @@
         @test promote(Stencil(1,1;center=1), Stencil(2.,2.;center=2)) == (Stencil(1.,1.;center=1), Stencil(2.,2.;center=2))
     end
 
+    @testset "apply_stencil" begin
+        v = [1, 2, 4, 8, 16, 32, 64, 128]
+        s = Stencil(1,2,3,4, center = 2)
+        @test apply_stencil(s,v, 2) == v[1] + 2*v[2] + 3*v[3] + 4*v[4]
+        @test apply_stencil(s,v, 4) == v[3] + 2*v[4] + 3*v[5] + 4*v[6]
+        @test apply_stencil_backwards(s,v, 3) == 4*v[1] + 3*v[2] + 2*v[3] + 1*v[4]
+        @test apply_stencil_backwards(s,v, 7) == 4*v[5] + 3*v[6] + 2*v[7] + 1*v[8]
+        @test apply_stencil(s,v, 2) isa Int
+        @test apply_stencil_backwards(s,v, 7) isa Int
+
+        v = [1, 2, 4, 8, 16, 32, 64, 128]
+        s = Stencil(1.,2.,3.,4., center = 2)
+        @test apply_stencil(s,v, 4) == v[3] + 2. *v[4] + 3. *v[5] + 4. *v[6]
+        @test apply_stencil_backwards(s,v, 7) == 4. *v[5] + 3. *v[6] + 2. *v[7] + v[8]
+        @test apply_stencil(s,v, 2) isa Float64
+        @test apply_stencil_backwards(s,v, 7) isa Float64
+
+        v = [1., 2., 4., 8., 16., 32., 64., 128.]
+        s = Stencil(1,2,3,4, center = 2)
+        @test apply_stencil(s,v, 2) == v[1] + 2*v[2] + 3*v[3] + 4*v[4]
+        @test apply_stencil_backwards(s,v, 3) == 4*v[1] + 3*v[2] + 2*v[3] + 1*v[4]
+        @test apply_stencil(s,v, 2) isa Float64
+        @test apply_stencil_backwards(s,v, 3) isa Float64
+
+        v = [@SVector[1, 2], @SVector[3, 4], @SVector[5, 6], @SVector[7, 8]]
+        s = Stencil(1,2, center = 1)
+        @test apply_stencil(s,v,1) == @SVector[7, 10]
+        @test apply_stencil_backwards(s,v,3) == @SVector[11, 14]
+        @test apply_stencil(s,v,1) isa SVector{2, Int}
+        @test apply_stencil_backwards(s,v,3) isa SVector{2, Int}
+
+        v = [@SVector[1., 2.], @SVector[3., 4.], @SVector[5., 6.], @SVector[7., 8.]]
+        s = Stencil(1,2, center = 1)
+        @test apply_stencil(s,v,1) == @SVector[7., 10.]
+        @test apply_stencil_backwards(s,v,3) == @SVector[11., 14.]
+        @test apply_stencil(s,v,1) isa SVector{2, Float64}
+        @test apply_stencil_backwards(s,v,3) isa SVector{2, Float64}
+
+        v = [@SVector[1, 2], @SVector[3, 4], @SVector[5, 6], @SVector[7, 8]]
+        s = Stencil(1.,2., center = 1)
+        @test apply_stencil(s,v,1) == @SVector[7., 10.]
+        @test apply_stencil_backwards(s,v,3) == @SVector[11., 14.]
+        @test apply_stencil(s,v,1) isa SVector{2, Float64}
+        @test apply_stencil_backwards(s,v,3) isa SVector{2, Float64}
+    end
+
     @testset "type stability" begin
         s_int = CenteredStencil(1,2,3)
         s_float = CenteredStencil(1.,2.,3.)
@@ -155,6 +203,27 @@
 
         @test SbpOperators.apply_stencil(ns, c, v, 4) == 5*7 + 11*11 + 6*13
         @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) == -3*3 - 7*5 - 4*7
+
+        # Different types in vector and stencil
+        ns = NestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-2.,2.), center=2)
+        @test SbpOperators.apply_inner_stencils(ns, c, 4) isa Stencil{Float64, 3}
+        @test SbpOperators.apply_inner_stencils(ns, c, 4) == Stencil(4.,9.,10.; center=2)
+        @test SbpOperators.apply_inner_stencils_backwards(ns, c, 4) isa Stencil{Float64, 3}
+        @test SbpOperators.apply_inner_stencils_backwards(ns, c, 4) == Stencil(-5.,-9.,-8.; center=2)
+
+        @test SbpOperators.apply_stencil(ns, c, v, 4) isa Float64
+        @test SbpOperators.apply_stencil(ns, c, v, 4) == 193.
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) isa Float64
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) == -158.
+
+        # Arrays of vectors
+        ns = NestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-2.,2.), center=2)
+        c = [  1,  3,  6, 10]
+        v = [@SVector[1, 2], @SVector[3, 4], @SVector[5, 6], @SVector[7, 8]]
+        @test SbpOperators.apply_stencil(ns, c, v, 2) isa SVector{2,Float64}
+        @test SbpOperators.apply_stencil(ns, c, v, 2) == 2v[1] + 5v[2] + 6v[3]
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 2) isa SVector{2,Float64}
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 2) == -4v[1] - 5v[2] - 3v[3]
     end
 
     @testset "type stability" begin