The following matrix equation $$B\Sigma + \Sigma B^T + C = 0$$ in $\Sigma$ $-$ for given $B$ and $C$ matrices $-$ appears in my work as a characterization of a covariance matrix. I have learned that this equation is known, in particularly in continuous time control theory, as Lyapunov's equation, and that there are various well known algorithms for solving it that exploit the special nature of this linear equation.
From googling I have also learned that there exist Matlab and Fortran implementations. I have found SLICOT and RECSY. Due to licensing issues access to SLICOT source has been stopped, though.
Most of my work is implemented in R, and as I have been unable to find an R interface to a solver, I consider writing one myself. My question is then if SLICOT is the best available Fortran (or C) library with an implementation of a solver of Lyapunov's equation? I am also interested in implementations that can handle large sparse $B$ matrices.