Merge pull request #15 from kalmarek/mk/reduce_allocations_arithmetic

reduce allocations in arithmetic

.github/workflows/CI.yml

@@ -1,15 +1,15 @@
 name: CI
 on:
   - pull_request
+  - push
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
     runs-on: ${{ matrix.os }}
     strategy:
-      fail-fast: false
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1'
           - 'nightly'
         os:
@@ -18,13 +18,16 @@ jobs:
           - windows-latest
         arch:
           - x64
+        allow_failures:
+          - julia: nightly
+      fail-fast: false
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v3
         env:
           cache-name: cache-artifacts
         with:
@@ -34,12 +37,15 @@ jobs:
             ${{ runner.os }}-test-${{ env.cache-name }}-
             ${{ runner.os }}-test-
             ${{ runner.os }}-
-      - uses: julia-actions/julia-buildpkg@v1
-      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-buildpkg@latest
+      - uses: julia-actions/julia-runtest@latest
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v2
         with:
-          file: lcov.info
+          file: ./lcov.info
+          flags: unittests
+          name: codecov-umbrella
+          fail_ci_if_error: false
   docs:
     name: Documentation
     runs-on: ubuntu-latest

.github/workflows/TagBot.yml

@@ -4,6 +4,11 @@ on:
     types:
       - created
   workflow_dispatch:
+    inputs:
+      lookback:
+        default: 3
+permissions:
+  contents: write
 jobs:
   TagBot:
     if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'

Project.toml

@@ -1,21 +1,22 @@
 name = "StarAlgebras"
 uuid = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 authors = ["Marek Kaluba <[email protected]>"]
-version = "0.1.7"
+version = "0.1.8"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-GroupsCore = "0.4"
-julia = "1"
+julia = "1.6"
 
 [extras]
-Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
 GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
+PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Pkg", "GroupsCore", "Random", "Test"]
+test = ["Pkg", "Groups", "GroupsCore", "PermutationGroups", "Random", "Test"]

src/algebra_elts.jl

@@ -50,3 +50,9 @@ aug(a::AlgebraElement) = sum(coeffs(a))
 LinearAlgebra.dot(a::AlgebraElement, v::AbstractVector) =
     LinearAlgebra.dot(StarAlgebras.coeffs(a), v)
 LinearAlgebra.dot(v::AbstractVector, a::AlgebraElement) = LinearAlgebra.dot(a, v)
+
+Base.copy(a::AlgebraElement) = AlgebraElement(copy(coeffs(a)), parent(a))
+function Base.deepcopy_internal(a::AlgebraElement, stackdict::IdDict)
+    haskey(stackdict, a) && return stackdict[a]
+    return AlgebraElement(deepcopy(coeffs(a)), parent(a))
+end

src/arithmetic.jl

@@ -6,6 +6,7 @@ Base.:(/)(X::AlgebraElement, a::Number) = inv(a) * X
 
 # TODO: handle this through mul!?
 Base.:(//)(X::AlgebraElement, a::Number) = AlgebraElement(coeffs(X) .// a, parent(X))
+Base.:div(X::AlgebraElement, a::Number) = AlgebraElement(div.(coeffs(X), a), parent(X))
 
 # ring structure:
 Base.:-(X::AlgebraElement) = neg!(similar(X), X)
@@ -33,28 +34,48 @@ end
 
 function neg!(res::AlgebraElement, X::AlgebraElement)
     @assert parent(res) === parent(X)
-    res.coeffs .= -coeffs(X)
+    res.coeffs .= .-coeffs(X)
     return res
 end
 
 function add!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement)
-    @assert parent(res) === parent(X) === parent(Y)
-    # res = (res === X || res === Y) ? similar(res) : res
-    res.coeffs .= coeffs(X) .+ coeffs(Y)
+    @assert parent(res) === parent(X)
+    @assert parent(X) === parent(Y)
+    if res === X
+        for (idx, y) in _nzpairs(coeffs(Y))
+            res[idx] += y
+        end
+    elseif res === Y
+        for (idx, x) in _nzpairs(coeffs(X))
+            res[idx] += x
+        end
+    else
+        zero!(res)
+        for (idx, x) in _nzpairs(coeffs(X))
+            res[idx] += x
+        end
+        for (idx, y) in _nzpairs(coeffs(Y))
+            res[idx] += y
+        end
+    end
     return res
 end
 
 function sub!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement)
     @assert parent(res) === parent(X) === parent(Y)
-    # res = (res === X || res === Y) ? similar(res) : res
-    res.coeffs .= coeffs(X) .- coeffs(Y)
+    neg!(res, Y)
+    add!(res, res, X)
     return res
 end
 
 function mul!(res::AlgebraElement, X::AlgebraElement, a::Number)
     @assert parent(res) === parent(X)
-    # res = (res === X) ? similar(res) : res
-    res.coeffs .= a .* coeffs(X)
+    if res !== X
+        zero!(res)
+    end
+    for (idx, x) in _nzpairs(coeffs(X))
+        res.coeffs[idx] = a * x
+    end
     return res
 end
 
@@ -69,12 +90,12 @@ function mul!(
 end
 
 function mul!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement)
-    @assert parent(res) === parent(X) === parent(Y)
     res = (res === X || res === Y) ? zero(res) : zero!(res)
     return fmac!(res, X, Y)
 end
 
 function fmac!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement)
+    @assert parent(res) === parent(X) === parent(Y)
     fmac!(coeffs(res), coeffs(X), coeffs(Y), parent(res).mstructure)
     return res
 end
@@ -91,9 +112,33 @@ function fmac!(
 )
     @assert res !== X
     @assert res !== Y
-    for (j, y) in _nzpairs(Y)
-        for (i, x) in _nzpairs(X)
-            res[mstr[i, j]] += x * y
+    lx, ly = size(mstr)
+    @assert all(iszero, @view(X[lx+1:end]))
+    @assert all(iszero, @view(Y[ly+1:end]))
+    for iy in 1:ly
+        y = Y[iy]
+        iszero(y) && continue
+        for ix in 1:lx
+            x = X[ix]
+            iszero(x) && continue
+            res[mstr[ix, iy]] += X[ix] * y
+        end
+    end
+    return res
+end
+
+function fmac!(
+    res::AbstractVector,
+    X::AbstractSparseVector,
+    Y::AbstractSparseVector,
+    mstr::MultiplicativeStructure,
+)
+    @assert res !== X
+    @assert res !== Y
+    for iy in SparseArrays.nonzeroinds(Y)
+        y = Y[iy]
+        for ix in SparseArrays.nonzeroinds(X)
+            res[mstr[ix, iy]] += X[ix] * y
         end
     end
     return res

src/mtables.jl

@@ -31,7 +31,7 @@ end
 
 MTable{Tw}(mt::AbstractMatrix{<:Integer}) where {Tw} = MTable{eltype(mt),Tw,typeof(mt)}(mt)
 
-MTable(b::AbstractBasis; table_size) = MTable{false}(b; table_size = table_size)
+MTable(b::AbstractBasis; table_size) = MTable{false}(b; table_size=table_size)
 
 function MTable{Tw}(basis::AbstractBasis; table_size) where {Tw}
     @assert length(table_size) == 2
@@ -78,7 +78,7 @@ struct CachedMTable{T,I,B<:Basis{T,I},M,Twisted} <: AbstractMTable{I,Twisted}
 end
 
 CachedMTable(basis::AbstractBasis; table_size) =
-    CachedMTable{false}(basis; table_size = table_size)
+    CachedMTable{false}(basis; table_size=table_size)
 
 function CachedMTable{Tw}(basis::AbstractBasis{T,I}; table_size) where {Tw,T,I}
     return CachedMTable{Tw}(basis, zeros(I, table_size))

src/show.jl

@@ -10,9 +10,9 @@ __needs_parens(a::AlgebraElement) = true
 
 function _coeff_elt_print(io, c, elt)
     print(io, c, "·")
-    __needs_parens(elt) && print(io, "(")
+    __needs_parens(elt) && print(io, '(')
     print(io, elt)
-    __needs_parens(elt) && print(io, ")")
+    __needs_parens(elt) && print(io, ')')
     return
 end
 
@@ -33,10 +33,12 @@ function Base.show(io::IO, a::AlgebraElement)
             if counter == 1
                 _coeff_elt_print(io, c, elt)
             else
-                print(io, ' ')
-                __prints_with_minus(c) || print(io, '+')
+                if __prints_with_minus(c)
+                    print(io, ' ')
+                else
+                    print(io, ' ', '+')
+                end
                 _coeff_elt_print(io, c, elt)
-                counter == length(nzeros) || print(io, ' ')
             end
         end
     else

src/types.jl

@@ -43,7 +43,7 @@ function StarAlgebra{Tw}(
     cache_size::Tuple{<:Integer,Integer};
     precompute=false
 ) where {Tw}
-    mstr = CachedMTable{Tw}(basis, table_size = cache_size)
+    mstr = CachedMTable{Tw}(basis, table_size=cache_size)
     precompute && complete!(mstr)
     return StarAlgebra(obj, basis, mstr)
 end
@@ -73,7 +73,7 @@ Base.eltype(a::AlgebraElement) = eltype(coeffs(a))
 
 ### constructing elements
 
-function Base.zero(A::AbstractStarAlgebra, T = Int)
+function Base.zero(A::AbstractStarAlgebra, T=Int)
     if hasbasis(A)
         I = SparseArrays.indtype(basis(A))
         return AlgebraElement(sparsevec(I[], T[], length(basis(A))), A)
@@ -83,7 +83,7 @@ function Base.zero(A::AbstractStarAlgebra, T = Int)
     )
 end
 
-function Base.one(A::AbstractStarAlgebra, T = Int)
+function Base.one(A::AbstractStarAlgebra, T=Int)
     hasbasis(A) && return A(one(object(A)), T)
     throw(
         "Algebra without basis; to construct one use the `AlgebraElement` constructor directly.",
@@ -101,31 +101,29 @@ function Base.isone(a::AlgebraElement)
     isone(a[k]) || return false
     return isone(b[k])
 end
-let SA = @static VERSION < v"1.3.0" ? :StarAlgebra : :AbstractStarAlgebra
-    @eval begin
-        function (A::$SA{O,T})(elt::T, S = Int) where {O,T}
-            if hasbasis(A)
-                b = basis(A)
-                i = b[elt]
-                return AlgebraElement(sparsevec([i], [one(S)], length(b)), A)
-            else
-                throw("Algebra without basis: cannot coerce $elt.")
-            end
+@eval begin
+    function (A::AbstractStarAlgebra{O,T})(elt::T, S=Int) where {O,T}
+        if hasbasis(A)
+            b = basis(A)
+            i = b[elt]
+            return AlgebraElement(sparsevec([i], [one(S)], length(b)), A)
+        else
+            throw("Algebra without basis: cannot coerce $elt.")
         end
+    end
 
-        function (A::$SA)(x::Number)
-            g = one(object(A))
-            res = A(g, typeof(x))
-            res[g] *= x
-            return res
-        end
+    function (A::AbstractStarAlgebra)(x::Number)
+        g = one(object(A))
+        res = A(g, typeof(x))
+        res[g] *= x
+        return res
     end
 end
 
-Base.similar(X::AlgebraElement, ::Type{T} = eltype(X)) where {T} =
+Base.similar(X::AlgebraElement, ::Type{T}=eltype(X)) where {T} =
     AlgebraElement(similar(coeffs(X), T), parent(X))
 
-function AlgebraElement{T}(X::AlgebraElement) where T
+function AlgebraElement{T}(X::AlgebraElement) where {T}
     v = coeffs(X)
     w = similar(v, T)
     w .= v