Solve matrix equation $X A Y + Y A^T X = B$ with diagonal $X$ and $Y$ Given square matrices $A$ and $B$, whereas $B$ is symmetric.   
How to solve matrix equation (unknowns are $X$ and $Y$)
$$ 
X A Y + Y A^T X = B
$$ 
with diagonal $X$ and $Y$ (if a solution exists)? And for which $A$ and $B$ is there a solution?
From inspecting the diagonal component, we get
$$
B_{ii} = X_{ii} A_{ii} Y_{ii} + Y_{ii} A_{ii} X_{ii}
$$
such that 
$$
X_{ii} Y_{ii} = \frac{1}{2} \frac{B_{ii}}{A_{ii}}.
$$
 A: This isn't a complete answer but I think it might be useful. First observe that for any matrix, $e_i^T M e_j = M_{ij}$ and for a diagonal matrix $e_i^T D = D_{ii} e_i^T$ and $D e_j = D_{jj} e_j$ where $e_i$ is the vector with $i$th element is $1$ and all other elements are $0$. So, we have
$$e_i^T X A Y e_j + e_i^T Y A^T X e_j = e_i^T B e_j $$
which simplifies to
$$E_{ij} : A_{ij} X_i Y_j + A_{ji} X_j Y_i = B_{ij}$$
where $X_i, Y_i$ are the diagonal elements (for easier typing). Note that $E_{ij} \equiv E_{ji}$, so instead of $n^2$ equations, we have $n(n+1)/2$ equations to solve. But we only have $2n$ unknowns, so for $n > 3$ the number of equations exceed the number of unknowns.
Regardless, we can construct the following matrix equation:
$$
\underbrace{\begin{bmatrix}
2 A_{11} & 0 & 0 & 0 & 0 & \dots\\
0 & A_{12} & A_{21} & 0 & 0 & \dots\\
0 & 0 & 0 & A_{13} & A_{31} & \dots\\
\vdots & \vdots & \vdots & \vdots & \vdots & \ddots
\end{bmatrix}}_{:=M}
\underbrace{\begin{bmatrix}
X_1 Y_1 \\ X_1 Y_2 \\ X_2 Y_1 \\ X_1 Y_3 \\ X_3 Y_1 \\ \vdots
\end{bmatrix}}_{:=x} =
\underbrace{\begin{bmatrix}
B_{11} \\ B_{12} \\ B_{13} \\ \vdots
\end{bmatrix}}_{:=b}
$$
So one necessary condition for a solution to exist is $b \in \operatorname{Im}M$. However, this is not enough. But you can try to find a solution by gradient descent methods by selecting an initial guess and iterating through a solution after this point.
Edit. Example for $n=3$.
$$
\begin{bmatrix}
2 A_{11} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & A_{12} & A_{21} & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & A_{13} & A_{31} & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 2 A_{22} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & A_{23} & A_{32} & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 A_{33}
\end{bmatrix}
\begin{bmatrix}
X_1 Y_1 \\ X_1 Y_2 \\ X_2 Y_1 \\ X_1 Y_3 \\ X_3 Y_1 \\ X_2 Y_2 \\ X_2 Y_3 \\  X_3 Y_2 \\ X_3 Y_3
\end{bmatrix} =
\begin{bmatrix}
B_{11} \\ B_{12} \\ B_{13} \\ B_{22} \\ B_{23} \\ B_{33}
\end{bmatrix}
$$
A: Let $A=[a_{i,j}],B=[b_{i,j}],X=[v_i],Y=[w_i]$. When a solution exists, necessarily, 
i) for every $i$, $a_{i,i}=0 \implies b_{i,i}=0$. 
ii) for every $i,j$, $a_{i,j}=a_{j,i}=0 \implies b_{i,j}=b_{j,i}=0$.
Assume that we consider generic matrices $A,B$ (for example, randomly choose them). Then, with probability $1$, the $(a_{i,j}),(b_{i,j})$ are non-zero and the conditions i),ii) are fulfilled. In particular, the $(v_i),(w_i)$ are non zero and we can put $v_1=1$ (if $(X,Y)$ is a solution, then $(aX,1/aY)$ too).
Thus one considers a system of $\dfrac{n(n+1)}{2}+1$ equations in the $2n$ unknowns $(a_{i,j}),(b_{i,j})$. As noted by obareey, the system is overdetermined for $n$ large enough.
Indeed, we can show that the equations are algebraically independent (when $A,B$ are generic); then (if we want that there is at least one solution) we must find $r_n=\dfrac{n(n+1)}{2}+1-2n=\dfrac{(n-1)(n-2)}{2}$ relations linking the $(a_{i,j}),(b_{i,j})$.
For example, we randomly choose $B$ and $n^2-r_n$ entries of the matrix $A$. Then one has $r_n$ supplementary unknowns.
For $n=2$, no problem. For $n=3$, $r_3=1$. For $n=4$, $r_4=3$.
That follows is an example for $n=4$. $u1,u2,u3$ are the $3$ supplementary unknowns.

One of the solutions for $u1,u2,u3$

The associated solution $(X,Y)$.

A: Let's assume that the diagonals are nonzero, i.e., $B_{ii} \neq 0$ and $A_{ii} \neq 0$, else simply $X_{ii} = Y_{ii} = 0$. More compact, we can write $x = diag(X)$ and $y = diag(Y)$
Then with
$$
x_{i} \circ y_{i} = \frac{1}{2} \frac{B_{ii}}{A_{ii}} = d_{i}. 
$$
where $\circ$ is the Hadamard (=element-wise) product.
We can reformulate to (similar to obareey's answer)
$$ 
A_{ij} x_i y_j + A_{ji} x_j y_i = B_{ij}
$$ 
and subsitute
$$
A_{ij} x_i \frac{d_j}{x_j} + A_{ji} x_j \frac{d_i}{x_i}  = B_{ij}
$$
By $z_{ij} = x_i / x_j$, we have
$$
A_{ij} d_j z_{ij} + A_{ji} d_i z_{ij}^{-1} = B_{ij}
$$
And in matrix form:
$$
A D \circ Z + D A^{T} \circ Z^{\circ -1} = B
$$
where $^{\circ -1}$ is an element-wise inverse. Thus
$$
A D \circ Z^{\circ 2} - B \circ Z + D A^{T}  = 0
$$
which are $n^2$ quadratic functions.
The remaining question is which of the two solutions of each of these $n^2$ equations satisfies $Z = x x^{-T}$?
