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(