Implement sum

src/BandedMatrices.jl

@@ -7,7 +7,7 @@ import Base: axes, axes1, getproperty, getindex, setindex!, *, +, -, ==, <, <=,
                >=, /, \, adjoint, transpose, showerror, convert, size, view,
                unsafe_indices, first, last, size, length, unsafe_length, step, to_indices,
                to_index, show, fill!, similar, copy, promote_rule, real, imag,
-               copyto!, Array
+               copyto!, Array, sum
 
 using Base.Broadcast: AbstractArrayStyle, DefaultArrayStyle, Broadcasted
 import Base.Broadcast: BroadcastStyle, broadcasted
@@ -99,4 +99,5 @@ end
 include("precompile.jl")
 
 
+
 end #module

src/banded/BandedMatrix.jl

@@ -1007,3 +1007,59 @@
     l,u = bandwidths(A)
     _BandedMatrix(reshape(resize!(vec(copy(bandeddata(A))), (l+u+1)*m), l+u+1, m), n, l,u)
 end
+
+function sum(A::BandedMatrix; dims=:)
+    if(dims == :)
+        l, u = bandwidths(A)
+        ret = zero(eltype(A))
+        if(l + u < 0)
+            return ret
+        end
+        n, m = size(A)
+        #Only get nonempty bands
+        lower, upper = min(n-1, l), min(m-1, u)
+        for i = -lower:upper
+            ret += sum(A[band(i)])
+        end
+        ret
+    elseif(dims > 2)
+        A
+    elseif(dims == 2)
+        l, u = bandwidths(A)
+        n, m = size(A)
+        ret = zeros(eltype(A), (n, 1))
+        if(l + u < 0)
+            return ret
+        end
+        lower, upper = min(n-1, l), min(m-1, u)
+        for i = -lower:upper
+            b = A[band(i)]
+            if(i <= 0)
+                ret[1-i:length(b)-i, 1] += b
+            else
+                ret[1:length(b), 1] += b
+            end
+        end
+        ret
+    elseif(dims == 1)
+        l, u = bandwidths(A)
+        n, m = size(A)
+        ret = zeros(eltype(A), (1, m))
+        if(l + u < 0)
+            return ret
+        end
+        lower, upper = min(n-1, l), min(m-1, u)
+        for i=-lower:upper
+            b = A[band(i)]
+            if(i <= 0)
+                ret[1, 1:length(b)] += b
+            else
+                ret[1, i+1:i+length(b)] += b
+            end
+        end
+        ret
+    else
+        throw(ArgumentError("dimension must be ≥ 1, got $dims"))
+    end
+end
+

test/runtests.jl

@@ -26,3 +26,4 @@ include("test_symbanded.jl")
 include("test_tribanded.jl")
 include("test_interface.jl")
 include("test_miscs.jl")
+include("test_sum.jl")

test/test_sum.jl

@@ -0,0 +1,15 @@
+module TestSum
+
+using Test, BandedMatrices, Random
+
+r = brand(Float64,rand(1:10_000),rand(1:10_000),rand(-20:100),rand(-20:100))
+matr = Matrix(r)
+@testset "sum" begin
+    @test sum(r) ≈ sum(matr) atol = 1e-10
+    @test sum(r; dims=2) ≈ sum(matr; dims=2) atol = 1e-10
+    @test sum(r; dims=1) ≈ sum(matr; dims=1) atol = 1e-10
+    @test sum(r; dims=3) == r
+    @test_throws ArgumentError sum(r; dims=0)
+end
+
+end