Use eigendecomposition-based Sylvester solver for symmetric case

Replace LAPACK trsyl!-based solver with a direct eigendecomposition
approach when both arguments are the same Hermitian matrix (as in polar
pullbacks). This avoids LAPACKException(1) for close eigenvalues.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: lkdvos

src/common/pullbacks.jl

@@ -11,5 +11,22 @@ function iszerotangent end
 iszerotangent(::Any) = false
 iszerotangent(::Nothing) = true
 
-# fallback
-_sylvester(A, B, C) = LinearAlgebra.sylvester(A, B, C)
+# Solve the Sylvester equation A*X + X*B + C = 0.
+# When A === B (same Hermitian PD matrix, as in polar pullbacks), use an
+# eigendecomposition-based solver to avoid LAPACK's trsyl! failing with
+# LAPACKException(1) for close eigenvalues.
+function _sylvester(A, B, C)
+    if A === B
+        return _sylvester_symm(A, C)
+    end
+    return LinearAlgebra.sylvester(A, B, C)
+end
+
+function _sylvester_symm(P, C)
+    D, Q = LinearAlgebra.eigen(LinearAlgebra.Hermitian(P))
+    Y = Q' * C * Q
+    @inbounds for j in axes(Y, 2), i in axes(Y, 1)
+        Y[i, j] = -Y[i, j] / (D[i] + D[j])
+    end
+    return Q * Y * Q'
+end