From de95d4a2b9b32cd82ee2e533d0c9d3defbcc60b6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Mon, 12 Aug 2024 12:38:34 +0200
Subject: [PATCH] MacaulayNullspace from MomentMatrix (#86)

---
 src/null.jl | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/src/null.jl b/src/null.jl
index a76faf8..280ee31 100644
--- a/src/null.jl
+++ b/src/null.jl
@@ -16,10 +16,22 @@ struct MacaulayNullspace{T,MT<:AbstractMatrix{T},BT<:SA.ExplicitBasis}
     basis::BT
     accuracy::T
 end
+
 function MacaulayNullspace(matrix::AbstractMatrix{T}, basis) where {T}
     return MacaulayNullspace(matrix, basis, Base.rtoldefault(T))
 end
 
+function MacaulayNullspace(
+    ν::MomentMatrix,
+    rank_check::RankCheck,
+    ldlt::LowRankLDLTAlgorithm = SVDLDLT(),
+)
+    M = value_matrix(ν)
+    chol = low_rank_ldlt(M, ldlt, rank_check)
+    @assert size(chol.L, 1) == LinearAlgebra.checksquare(M)
+    return MacaulayNullspace(chol.L, ν.basis, accuracy(chol))
+end
+
 function Base.getindex(
     null::MacaulayNullspace{T,MT,<:MB.SubBasis{B}},
     monos,
@@ -69,8 +81,8 @@ function compute_support!(
     rank_check::RankCheck,
     solver::ImageSpaceSolver,
 )
-    ν.support =
-        solve(MacaulayNullspace(ν, rank_check, solver.ldlt), solver.null)
+    null = MacaulayNullspace(ν, rank_check, solver.ldlt)
+    ν.support = solve(null, solver.null)
     return
 end