LieAlgebras: reduce allocations (#4165)

.gitignore

@@ -4,6 +4,7 @@
 /docs/src/Nemo/
 /docs/src/Experimental/
 profile.pb.gz
+alloc-profile.pb.gz
 LocalPreferences.toml
 Manifest.toml
 .DS_Store

experimental/LieAlgebras/src/AbstractLieAlgebra.jl

@@ -265,9 +265,10 @@ function lie_algebra(
 
   npos = n_positive_roots(rs)
   # set Chevalley basis
-  r_plus = basis(L)[1:npos]
-  r_minus = basis(L)[(npos + 1):(2 * npos)]
-  h = basis(L)[(2 * npos + 1):end]
+  basis_L = basis(L)
+  r_plus = basis_L[1:npos]
+  r_minus = basis_L[(npos + 1):(2 * npos)]
+  h = basis_L[(2 * npos + 1):end]
   chev = (r_plus, r_minus, h)
   set_root_system_and_chevalley_basis!(L, rs, chev)
   return L
@@ -286,6 +287,7 @@ function _struct_consts(R::Field, rs::RootSystem, extraspecial_pair_signs)
   N = _N_matrix(rs, extraspecial_pair_signs)
 
   struct_consts = Matrix{sparse_row_type(R)}(undef, n, n)
+  beta_i_plus_beta_j = zero(RootSpaceElem, rs)
   for i in 1:nroots, j in i:nroots
     if i == j
       # [e_βi, e_βi] = 0
@@ -294,12 +296,14 @@ function _struct_consts(R::Field, rs::RootSystem, extraspecial_pair_signs)
     end
     beta_i = root(rs, i)
     beta_j = root(rs, j)
-    if iszero(beta_i + beta_j)
+    beta_i_plus_beta_j = add!(beta_i_plus_beta_j, beta_i, beta_j)
+    if iszero(beta_i_plus_beta_j)
       # [e_βi, e_-βi] = h_βi
       struct_consts[i, j] = sparse_row(
-        R, collect((nroots + 1):(nroots + nsimp)), _vec(coefficients(coroot(rs, i)))
+        R, collect((nroots + 1):(nroots + nsimp)), _vec(coefficients(coroot(rs, i)));
+        sort=false,
       )
-    elseif ((fl, k) = is_root_with_index(beta_i + beta_j); fl)
+    elseif ((fl, k) = is_root_with_index(beta_i_plus_beta_j); fl)
       # complicated case
       if i <= npos
         struct_consts[i, j] = sparse_row(R, [k], [N[i, j]])
@@ -313,15 +317,16 @@ function _struct_consts(R::Field, rs::RootSystem, extraspecial_pair_signs)
       struct_consts[i, j] = sparse_row(R)
     end
 
-    # # [e_βj, e_βi] = -[e_βi, e_βj]
+    # [e_βj, e_βi] = -[e_βi, e_βj]
     struct_consts[j, i] = -struct_consts[i, j]
   end
   for i in 1:nsimp, j in 1:nroots
     # [h_i, e_βj] = <β_j, α_i> e_βj
     struct_consts[nroots + i, j] = sparse_row(
       R,
       [j],
-      [dot(coefficients(root(rs, j)), view(cm, i, :))],
+      # [dot(coefficients(root(rs, j)), view(cm, i, :))], # currently the below is faster
+      [only(coefficients(root(rs, j)) * cm[i, :])],
     )
     # [e_βj, h_i] = -[h_i, e_βj]
     struct_consts[j, nroots + i] = -struct_consts[nroots + i, j]
@@ -352,25 +357,24 @@ function _N_matrix(rs::RootSystem, extraspecial_pair_signs::Vector{Bool})
         j -> !is_positive_root(alpha_i + beta_l - simple_root(rs, j)),
         1:(i - 1),
       ) || continue
-      p = 0
-      while is_root(beta_l - p * alpha_i)
-        p += 1
-      end
-      N[i, l] = (extraspecial_pair_signs[k - nsimp] ? 1 : -1) * p
+      p = _root_string_length_down(beta_l, alpha_i) + 1
+      N[i, l] = (extraspecial_pair_signs[k - nsimp] ? p : -p)
       N[l, i] = -N[i, l]
     end
   end
 
   # special pairs
+  alpha_i_plus_beta_j = zero(RootSpaceElem, rs)
   for (i, alpha_i) in enumerate(positive_roots(rs))
     for (j, beta_j) in enumerate(positive_roots(rs))
       i < j || continue
-      is_positive_root(alpha_i + beta_j) || continue
+      alpha_i_plus_beta_j = add!(alpha_i_plus_beta_j, alpha_i, beta_j)
+      is_positive_root(alpha_i_plus_beta_j) || continue
       l = findfirst(
-        l -> is_positive_root(alpha_i + beta_j - simple_root(rs, l)), 1:nsimp
+        l -> is_positive_root(alpha_i_plus_beta_j - simple_root(rs, l)), 1:nsimp
       )::Int
       l == i && continue # already extraspecial
-      fl, l_comp = is_positive_root_with_index(alpha_i + beta_j - simple_root(rs, l))
+      fl, l_comp = is_positive_root_with_index(alpha_i_plus_beta_j - simple_root(rs, l))
       @assert fl
       t1 = 0
       t2 = 0
@@ -383,10 +387,7 @@ function _N_matrix(rs::RootSystem, extraspecial_pair_signs::Vector{Bool})
         t2 = N[l, m] * N[j, m] * dot(root_m, root_m)//dot(alpha_i, alpha_i)
       end
       @assert t1 - t2 != 0
-      p = 0
-      while is_root(beta_j - p * alpha_i)
-        p += 1
-      end
+      p = _root_string_length_down(beta_j, alpha_i) + 1
       N[i, j] = Int(sign(t1 - t2) * sign(N[l, l_comp]) * p) # typo in CMT04
       N[j, i] = -N[i, j]
     end

experimental/LieAlgebras/src/LieAlgebra.jl

@@ -49,7 +49,9 @@ Return the `i`-th basis element of the Lie algebra `L`.
 function basis(L::LieAlgebra, i::Int)
   @req 1 <= i <= dim(L) "Index out of bounds."
   R = coefficient_ring(L)
-  return L([(j == i ? one(R) : zero(R)) for j in 1:dim(L)])
+  v = zero_matrix(R, 1, dim(L))
+  v[1, i] = one(R)
+  return L(v)
 end
 
 @doc raw"""
@@ -453,10 +455,12 @@ end
 
 function adjoint_matrix(x::LieAlgebraElem{C}) where {C<:FieldElem}
   L = parent(x)
-  return sum(
-    c * g for (c, g) in zip(coefficients(x), adjoint_matrices(L));
-    init=zero_matrix(coefficient_ring(L), dim(L), dim(L)),
-  )
+  mat = zero_matrix(coefficient_ring(L), dim(L), dim(L))
+  tmp = zero(mat)
+  for (c, g) in zip(coefficients(x), adjoint_matrices(L))
+    mat = addmul!(mat, g, c, tmp)
+  end
+  return mat
 end
 
 function _adjoint_matrix(S::LieSubalgebra{C}, x::LieAlgebraElem{C}) where {C<:FieldElem}

experimental/LieAlgebras/src/LieAlgebraIdeal.jl

@@ -270,7 +270,7 @@ where `L` is the Lie algebra where `I` lives in.
 function lie_algebra(I::LieAlgebraIdeal)
   LI = lie_algebra(basis(I))
   L = base_lie_algebra(I)
-  emb = hom(LI, L, basis(I))
+  emb = hom(LI, L, basis(I); check=false)
   return LI, emb
 end
 

experimental/LieAlgebras/src/LieAlgebraModule.jl

@@ -53,7 +53,9 @@ Return the `i`-th basis element of the Lie algebra module `V`.
 function basis(V::LieAlgebraModule, i::Int)
   @req 1 <= i <= dim(V) "Index out of bounds."
   R = coefficient_ring(V)
-  return V([(j == i ? one(R) : zero(R)) for j in 1:dim(V)])
+  v = zero_matrix(R, 1, dim(V))
+  v[1, i] = one(R)
+  return V(v)
 end
 
 @doc raw"""

experimental/LieAlgebras/src/LieAlgebras.jl

@@ -84,6 +84,7 @@ import ..Oscar:
   tensor_product,
   weyl_vector,
   word,
+  zero!,
   zero_map,
   ⊕,
   ⊗

experimental/LieAlgebras/src/LieSubalgebra.jl

@@ -301,7 +301,7 @@ where `L` is the Lie algebra where `S` lives in.
 function lie_algebra(S::LieSubalgebra)
   LS = lie_algebra(basis(S))
   L = base_lie_algebra(S)
-  emb = hom(LS, L, basis(S))
+  emb = hom(LS, L, basis(S); check=false)
   return LS, emb
 end
 

experimental/LieAlgebras/src/LinearLieAlgebra.jl

@@ -21,7 +21,7 @@ Return the basis `basis(L)` of the Lie algebra `L` in the underlying matrix
 representation.
 """
 function matrix_repr_basis(L::LinearLieAlgebra{C}) where {C<:FieldElem}
-  return Vector{dense_matrix_type(C)}(L.basis)
+  return L.basis::Vector{dense_matrix_type(C)}
 end
 
 @doc raw"""
@@ -31,7 +31,7 @@ Return the `i`-th element of the basis `basis(L)` of the Lie algebra `L` in the
 underlying matrix representation.
 """
 function matrix_repr_basis(L::LinearLieAlgebra{C}, i::Int) where {C<:FieldElem}
-  return (L.basis[i])::dense_matrix_type(C)
+  return matrix_repr_basis(L)[i]
 end
 
 ###############################################################################
@@ -130,10 +130,12 @@ Return the Lie algebra element `x` in the underlying matrix representation.
 """
 function Generic.matrix_repr(x::LinearLieAlgebraElem)
   L = parent(x)
-  return sum(
-    c * b for (c, b) in zip(_matrix(x), matrix_repr_basis(L));
-    init=zero_matrix(coefficient_ring(L), L.n, L.n),
-  )
+  mat = zero_matrix(coefficient_ring(L), L.n, L.n)
+  tmp = zero(mat)
+  for (c, b) in zip(coefficients(x), matrix_repr_basis(L))
+    mat = addmul!(mat, b, c, tmp)
+  end
+  return mat
 end
 
 function bracket(