Fix sampled matrix integration over dimensional axes

Project.toml

@@ -71,6 +71,7 @@ SciMLLogging = "1.10.1, 2"
 SciMLTesting = "1"
 StaticArrays = "1.9.8"
 Test = "1.10"
+Unitful = "1"
 Zygote = "0.7.10"
 julia = "1.10"
 
@@ -90,6 +91,7 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["ADTypes", "Arblib", "StaticArrays", "SafeTestsets", "SciMLTesting", "Test", "Distributions", "ChainRulesCore", "FastGaussQuadrature", "FastTanhSinhQuadrature", "Cuba", "Cubature", "MCIntegration", "DifferentiationInterface", "HAdaptiveIntegration"]
+test = ["ADTypes", "Arblib", "StaticArrays", "SafeTestsets", "SciMLTesting", "Test", "Distributions", "ChainRulesCore", "FastGaussQuadrature", "FastTanhSinhQuadrature", "Cuba", "Cubature", "MCIntegration", "DifferentiationInterface", "HAdaptiveIntegration", "Unitful"]

src/sampled.jl

@@ -72,8 +72,13 @@ function evalrule(data::AbstractMatrix{T}, weights, dim) where {T}
         n == 0 && throw(ArgumentError("No points to integrate"))
         if dim == 2 || dim == ndims(data)
             # Integration over columns (last dimension) - the common case
-            out = zeros(T, m)
-            @inbounds for i in 1:n
+            w1 = weights[1]
+            out_type = m == 0 ? typeof(w1 * zero(T)) : typeof(w1 * data[1, 1])
+            out = Vector{out_type}(undef, m)
+            @inbounds @simd for j in 1:m
+                out[j] = w1 * data[j, 1]
+            end
+            @inbounds for i in 2:n
                 w = weights[i]
                 @simd for j in 1:m
                     out[j] += w * data[j, i]

test/sampled_tests.jl

@@ -1,4 +1,4 @@
-using Integrals, Test
+using Integrals, Test, Unitful
 @testset "Sampled Integration" begin
     lb = 0.4
     ub = 1.1
@@ -32,6 +32,16 @@ using Integrals, Test
     end
 end
 
+@testset "Sampled integration with dimensional axes" begin
+    data = ones(2, 3)
+    axis = [0.0, 1.0, 2.0]u"m"
+
+    prob = SampledIntegralProblem(data, axis; dim = 2)
+    sol = solve(prob, TrapezoidalRule())
+
+    @test sol.u == [2.0, 2.0]u"m"
+end
+
 @testset "Caching interface" begin
     x = 0.0:0.1:1.0
     y = sin.(x)