From 1861eb7075a6963a04293934d209d0a13b915ee3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Wed, 27 Sep 2023 15:42:28 +0200
Subject: [PATCH] Fix edge case for block diagonalization

---
 src/Certificate/Symmetry/block_diag.jl | 18 ++++++++++--------
 src/Certificate/Symmetry/wedderburn.jl | 10 +++++-----
 test/symmetry.jl                       | 23 +++++++++++++++++++++++
 3 files changed, 38 insertions(+), 13 deletions(-)

diff --git a/src/Certificate/Symmetry/block_diag.jl b/src/Certificate/Symmetry/block_diag.jl
index 6b01d29ae..091b9477c 100644
--- a/src/Certificate/Symmetry/block_diag.jl
+++ b/src/Certificate/Symmetry/block_diag.jl
@@ -180,7 +180,7 @@ end
 
 function ordered_block_diag(As, d)
     U = block_diag(As, d)
-    U === nothing && return nothing
+    isnothing(U) && return nothing
     iU = U'
     @assert iU ≈ inv(U)
     Bs = [iU * A * U for A in As]
@@ -199,7 +199,8 @@ function ordered_block_diag(As, d)
         # should work, this trick is similar to [CGT97].
         #
         # [CGT97] Corless, R. M.; Gianni, P. M. & Trager, B. M.
-        # A reordered Schur factorization method for zero-dimensional polynomial systems with multiple roots Proceedings of the 1997 international symposium on Symbolic and algebraic computation,
+        # A reordered Schur factorization method for zero-dimensional polynomial systems with multiple roots
+        # Proceedings of the 1997 international symposium on Symbolic and algebraic computation,
         # 1997, 133-140
         R = sum(λ .* refs)
         C = sum(λ .* Cs)
@@ -247,13 +248,14 @@ function block_diag(As, d)
             offset = 0
             for v in values(blocks)
                 @assert iszero(length(v) % d)
-                if length(v) == d
-                    V = Z[:, v]
-                else
+                V = Z[:, v]
+                if length(v) != d
                     Cs = [B[v, v] for B in Bs]
-                    V = block_diag(Cs, d)
-                    V === nothing && break
-                    V *= transpose(Z[:, v])
+                    _V = block_diag(Cs, d)
+                    if isnothing(_V)
+                        break
+                    end
+                    V *= _V
                 end
                 U[:, offset.+eachindex(v)] = V
                 offset += length(v)
diff --git a/src/Certificate/Symmetry/wedderburn.jl b/src/Certificate/Symmetry/wedderburn.jl
index efe8f6cbb..71adffce5 100644
--- a/src/Certificate/Symmetry/wedderburn.jl
+++ b/src/Certificate/Symmetry/wedderburn.jl
@@ -185,14 +185,14 @@ function _gram_basis(pattern::Pattern, basis, ::Type{T}) where {T}
         if d > 1
             if m > 1
                 U = ordered_block_diag(S, d)
+                if isnothing(U)
+                    error(
+                        "Could not simultaneously block-diagonalize into $m identical $dx$d blocks",
+                    )
+                end
             else
                 U = Matrix{T}(LinearAlgebra.I, N, N)
             end
-            if U === nothing
-                error(
-                    "Could not simultaneously block-diagonalize into $m identical $dx$d blocks",
-                )
-            end
             # From Example 1.7.3 of
             # Sagan, The symmetric group, Springer Science & Business Media, 2001
             # we know that there exists `C` such that `Q = kron(C, I)` if we use
diff --git a/test/symmetry.jl b/test/symmetry.jl
index 1eca7f667..c267d8524 100644
--- a/test/symmetry.jl
+++ b/test/symmetry.jl
@@ -184,6 +184,29 @@ function test_block_diag()
     d = 2
     U = SumOfSquares.Certificate.Symmetry.ordered_block_diag(A, d)
     @test SumOfSquares.Certificate.Symmetry.ordered_block_check(U, A, d)
+    α = 0.75
+    A = [
+        Matrix{Int}(I, 6, 6),
+        [
+            1 0 0 0 0 0
+            0 -1 0 0 0 0
+            0 0 0 0 1 0
+            0 0 0 0 0 -1
+            0 0 1 0 0 0
+            0 0 0 -1 0 0
+        ],
+        [
+            -0.5 α 0 0 0 0
+            -α -0.5 0 0 0 0
+            0 0 -0.5 α 0 0
+            0 0 -α -0.5 0 0
+            0 0 0 0 -0.5 α
+            0 0 0 0 -α -0.5
+        ],
+    ]
+    d = 2
+    U = SumOfSquares.Certificate.Symmetry.ordered_block_diag(A, d)
+    @test SumOfSquares.Certificate.Symmetry.ordered_block_check(U, A, d)
 end
 
 function runtests()