I am trying to find the symmetric solution $X\in \mathbb{R}^{p\times p}$ of following matrix equation:

$AXB + (AXB)^T + cX = D$

where $A,B\in \mathbb{R}^{p\times p}$ are symmetric positive definite matrix, and $D\in \mathbb{R}^{p\times p}$ will also be symmetric.

e.g.

if $c=0$ and one of $A$ or $B$ equals identity matrix $I$, we call this equation Sylvester Equation, and have closed form solution using SVD.

The existence and uniqueness of the solution can be proved under the condition that $\beta_i\alpha_j + \alpha_i\beta_j + c \neq 0$ for all $i,j=1,\cdots ,p$ where $\alpha_i$ and $\beta_i$ are eigenvalue of $A,B$ respectively.

Currently, I can't find closed form solution of this equation when $A,B,D$ are arbitrary matrix. It seems like the only way to solve it is transforming this equation into linear equation using kronecker product

$(B\otimes A + A\otimes B + cI)Vec(X) = Vec(D)$

and solve this problem using conjugate gradient descent naively. but it's computational unafforfable when the size $p$ is relative large, since we need to solve a $p^2 \times p^2$ problem.

What if we take advantage of the special structure of the matrix? e.g. $X$ is a symmetric matrix and $A, B$ are symmetric positive definite matrix.

Can someone figure it out? Thanks!