Consider the 6th-order tensor ${\mathbb M}$, whose components ${\mathbb M}_{ijklmn}$ are unity if $\,(i=k=m)$ and $(j=l=n),\,$ but zero otherwise.
This tensor can be used, along with the Frobenius Product (:), to rewrite the Hadamard Product ($\circ$) between two matrices as
$$\eqalign{
B\circ A &= B:{\mathbb M}:A \cr\cr
}$$
This makes the finding the gradient (a 4th-order tensor) especially simple
$$\eqalign{
F &= A\circ B = B\circ A = B:{\mathbb M}:A \cr\cr
dF &= B:{\mathbb M}:dA \cr\cr
\frac{\partial F}{\partial A} &= B:{\mathbb M} \cr
}$$
According to the comments the only ugly part was the devectorization.
Consider the 3x3 matrix $${\bf X}=\left[\begin{array}{ccc}
1&2&3\\4&5&6\\7&8&9\end{array}\right]$$
It's lexical vectorization is:
$$\text{vec}({\bf X}) = \left[\begin{array}{ccccccccc}1&4&7&2&5&8&3&6&9\end{array}\right]^T$$
We consider that kind of vectorization here.
The Matlab / Octave expression:
kron([1,0,0]',eye(3))'*kron(R(:),[1,0,0])+
kron([0,1,0]',eye(3))'*kron(R(:),[0,1,0])+
kron([0,0,1]',eye(3))'*kron(R(:),[0,0,1])-R
evaluates to the zero matrix for several random matrices R, each new term "selecting" a new row.
Then a guess would be that
$${\bf R}=\sum_k({\bf v_k} \otimes {\bf I})^T(\text{vec}({\bf R})\otimes {\bf v_k}^T)$$
With $\bf v_k$ being the Dirac or selector vector:
$$({\bf v_k})_i = \cases{0, k \neq i\\1, k=i}$$
Feel free to try and simplify it if you want.
Best Answer
We define a function as $f: \mathbb{R}^k \to \mathbb{R}^{k \times k}, f(x)=(xx^T)\circ A$.
The domain of the function are vectors in $\mathbb{R}^k$ and the range of the function are matrices in $\mathbb{R}^{k \times k}$.
We want to evaluate the derivative of the function with respect to $x$. For each entries of the matrix, we can calculate the partial derivative of that entry with respect to the the different entries in the vector $x$. The result can therefore be organized into a 3rd order tensor of dimension $k \times k \times k$ that we can denote by $\frac{\partial f(x)}{\partial x}$ where the partial derivative of the $i,j$ entry of $f(x)$ with respect to $k^{th}$ entry in $x$ is denoted by $\frac{\partial f(x)}{\partial x}_{i,j,k}$. We know that the $i,j$ entry of $f(x)_{i,j}=x_ix_ja_{ij}$, so $$\frac{\partial f(x)}{\partial x}_{ijk}=\begin{cases} a_{ij}x_j, s=i\\ a_{ij}x_i, s=j\\ 0 , \text{ otherwise} \end{cases}$$