Implementing the Bartels–Stewart algorithm in Eigen3 -- real matrices only?

Based off this question and solution -- Implementing the Bartels–Stewart algorithm in Eigen3? -- I am trying to solve Lyapunov equations (AX + XA^T = C) using the Eigen library, but am limited to real matrices.

The R (with c++) code below works, but involves complex numbers. It can definitely be simplified (since in this framing, there is no B matrix), but the main difficulty is the reliance on complex numbers. The real schur form seems to be the standard alternative in this case, but the Eigen function matrix_function_solve_triangular_sylvester then does not work because the input matrix is not upper triangular, but is upper block triangular. I would be happy to see suggestions to a) remove the need for complex numbers, and then if that is possible, b) any efficiency improvements.

library(expm)
library(Rcpp)
library(RcppEigen)
library(inline)
# R -----------------------------------------------------------------------
d<-6 #dimensions 
A<-matrix(rnorm(d^2),d,d) #continuous time transition
G <- matrix(rnorm(d^2),d,d)
C<-G %*% t(G) #continuous time pos def error
AHATCH<-A %x% diag(d) + diag(d) %x% A
Xtrue<-matrix(-solve(AHATCH,c(C)), d) #asymptotic error from continuous time


# c++ in R ---------------------------------------------------------------------

sylcpp <- '

    using Eigen::Map;
    using Eigen::MatrixXd;

  // Map the double matrix A from Ar
   const Map<MatrixXd> A(as<Map<MatrixXd> >(Ar));

   // Map the double matrix Q from Qr
   const Map<MatrixXd> Q(as<Map<MatrixXd> >(Qr));


  Eigen::MatrixXd B = A.transpose();

  Eigen::ComplexSchur<Eigen::MatrixXd> SchurA(A);
  Eigen::MatrixXcd R = SchurA.matrixT();
  Eigen::MatrixXcd U = SchurA.matrixU();

  Eigen::ComplexSchur<Eigen::MatrixXd> SchurB(B);
  Eigen::MatrixXcd S = SchurB.matrixT();
  Eigen::MatrixXcd V = SchurB.matrixU();

  Eigen::MatrixXcd F = (U.adjoint() * Q) * V;

  Eigen::MatrixXcd Y = Eigen::internal::matrix_function_solve_triangular_sylvester(R, S, F);

  Eigen::MatrixXd X = ((U * Y) * V.adjoint()).real();

  return wrap(X);
'



syl <- cxxfunction(signature(Ar = "matrix",Qr='matrix'), sylcpp, plugin = "RcppEigen")

X=syl(A,-C)
X-Xtrue #approx zero