After hustling with this identity for a bit, this is what I was able to come up with. First thing to pay attention to is that $\nabla \cdot (\vec A\times \vec B )$ is the divergence of the cross product vector field.
The interpretation for the cross product vector field depends on the domain of the problem, but we can abstract that and just think that $\vec A\times \vec B$ generates a new vector field $\vec C$.
Hence, $\nabla \cdot \vec C$ can be understood using the common intuition for divergence, which is that it measures the amount of "fluid" being removed or added to the field (sink/source). For a better explanation on this, you can check out this videos by ThreeBrownOneBlue.
Now, to get a feel for what is going on in the identity, we will explore a very simple 3D case in which the vector field $\vec B$ varies only in the $\hat k$ axis, and $\vec A$ is equal to 1 in the $\hat i$ axis.
$$\vec A(x,y,z)= 1\hat i + 0 \hat j + 0\hat k$$
$$\vec B(x,y,z)= 0\hat i + (z+2) \hat j + 0\hat k$$
Below we draw a cube of 1x1x1, and show in a diagram the how the vector fields acts.
$$\\$$
Note that $\nabla \cdot \vec C$ is the difference between the flow leaving the top face and the flow entering the bottom face. From this we can start to understand how the curl of $\vec B$ affects the divergence of $\vec C$. The $\nabla \times \vec B$ is commonly viewed as the "rotation" of the field, but, such rotation implies that field $\vec B$ is changing in the perpendicular direction, which is exactly where $\vec A \times \vec B$ will be pointing. Hence, a curl in $\vec B$ means an acceleration in $\vec A \times \vec B$ and therefore, a variation in the divergence.
In this example, we can see that the divergence in the cross product field can be calculated in two ways. Either evaluating the change going inside one face and leaving other. Or by calculating the variation occurring inside the cube. This two ways of calculating are the equivalence shown in the identity.
Therefore, calculating the divergence using the surfaces we get
$$\nabla\cdot \vec C = +1\cdot A - 2 \cdot A = -1 \cdot A = -1$$
Where $A$ is the area of the face cube. And using the curl, we get;
$$\nabla\cdot \vec C = \vec B \cdot (\nabla \times \vec A) \cdot h - \vec A \cdot (\nabla \times \vec B) = 0 - 1\hat i \cdot (1\hat i) \cdot h = -1$$
Where $h$ is the height of the cube.
We calculated variation in the cube with 1x1x1, but the same could've been done for the infinitesimally small cube, which would eliminate $A$ and $h$ from the calculation.
Finally, the case illustrated is for very simple fields, but more complex fields can be thought in the same way, by just breaking the effect in each direction and applying the same logic. I also hope that it is clear that the same explanation shows how the curl of $\vec A$ appears in the equation. If $\vec A$ had the same curl as $\vec B$ $$\vec A(x,y,z) = (z+2)\hat i + 0 \hat j + 0 \hat k$$
Then, the curl of $\vec A$ would also "accelerate" the cross product field in the $\hat k$ axis, decreasing the divergence even more.
Hopefully this gives you some intuition on the identity. If something is not clear let me know, and I can try to do some more diagrams.
Due to the complexity of the topic of the question, it is not possible to expound here all the associated analytic developments, therefore I try to answer by giving precise references and surveying the results presented there.
- Do solutions for $f$ exist?
Yes: to see this, note that, after formally manipulating your equation (perhaps going backward in the steps of your deduction), we get
$$
\begin{split}
\nabla\cdot\nabla f(\vec{r})&=\left( \frac{1}{2}\vec{B}\times\vec{r} - \nabla f(\vec{r}) \right)\cdot\frac{\nabla \rho(\vec{r})}{\rho(\vec{r})}\\
\rho(\vec{r})\nabla\cdot\nabla f(\vec{r})&=\left( \frac{1}{2}\vec{B}\times\vec{r} - \nabla f(\vec{r}) \right)\cdot\nabla \rho(\vec{r})\\
\rho(\vec{r})\nabla\cdot\nabla f(\vec{r})&=\left( \frac{1}{2}\vec{B}\times\vec{r} - \nabla f(\vec{r}) \right)\cdot\nabla \rho(\vec{r})\\
\rho(\vec{r})\nabla\cdot\nabla f(\vec{r})+\nabla \rho(\vec{r})\cdot\nabla f(\vec{r})&=\frac{1}{2}\vec{B}\times\vec{r}\cdot\nabla \rho(\vec{r})\\
\end{split}
$$
i.e.
$$
\nabla\cdot\big(\rho(\vec{r})\nabla f(\vec{r})\big)=\frac{1}{2}\vec{B}\times\vec{r}\cdot\nabla \rho(\vec{r}).\label{1}\tag{1}
$$
Equation \eqref{1} is a standard linear elliptic equations in divergence form, and as such it admits a fundamental solution (and also a Green function) i.e. a solution of the equation
$$
\nabla\cdot\big(\rho(\vec{r})\nabla \mathscr{E}(\vec{r},\vec{s})\big)=\delta(\vec{r}-\vec{s})\label{2}\tag{2}
$$
where $\delta$ is the usual Dirac distribution, under fairly general conditions. In particular $\rho(\vec{r})$ is required only to be a bounded and measurable function: moreover, the problem has a solution even if, instead of a scalar function $\rho(\vec{r})$, we have a second-order tensor function
$$
\vec{r}\mapsto\overline{\overline{\boldsymbol{\rho}}}(\vec{r})\in\Bbb R^{3\times 3}\quad\vec{r}\in\Bbb R^3 \label{3}\tag{3}
$$ with bounded measurable components. This result is due to Walter Littman, Guido Stampacchia and Hans Weinberger: the paper [2] (where the original result is to be found) is not an easy read and even the used notation is not a modern one, while the set of lecture notes [3] is a better readable source (with also an updated notation) but is written in French. A more modern reference is [1], where the authors prove also that the tensor of the coefficients \eqref{3} can be asymmetric.
- Is it possible to say something else than the above equation on the form of the solutions $f$? (Especially: is there an analytical solution?)
Yes (with some cautions): for the general case, several properties of the solution to equation \eqref{2}, and these properties "interact" with the ones of the given function $\frac{1}{2}\vec{B}\times\vec{r}\cdot\nabla\rho(\vec{r})$ (have again a look at the items in the bibliography for the details and also at this answer). For example it is known that the quotient between $\mathscr{E}(\vec{r},\vec{s})$ and the fundamental solution of Laplace equation is bounded by two positive constants $K$ and $K^{-1}$. And also you can use $\mathscr{E}(\vec{r},\vec{s})$ to construct a fairly explicit general solution of \eqref{1} as a convolution integral:
$$
f(\vec{r})=\int\limits_{\Bbb R^3}\mathscr{E}(\vec{r},\vec{s})\big(\vec{B}\times\vec{s}\cdot\nabla\rho(\vec{s})\big)\mathrm{d}\vec{s}
$$
However, it is not always possible to construct explicitly the fundamental solution $\mathscr{E}(\vec{r},\vec{s})$ (in this respect see also this answer) thus, apart from the case $\rho(\vec{r})\equiv\mathrm{const}.$ where \eqref{1} reduce to the Laplace equation
$$
\nabla\cdot\nabla f(\vec{r})=0,
$$
I am not aware of the existence of any exact (meant as expressed by means of more or less elementary functions) solutions to \eqref{1}.
Bibliography
[1] Michael Grüter and Kjell-Ove Widman (1982), "The Green function for uniformly elliptic equations", Manuscripta Mathematica 37, pp. 303-342 Zbl 0485.35031.
[2] Walter Littman, Hans Weinberger and Guido Stampacchia (1962), "Regular points for elliptic equations with discontinuous coefficients", Annali della Scuola Normale Superiore di Pisa, Classe di Scienze, serie III, Vol. 17, n° 1-2, pp. 43-77, MR161019, Zbl0116.30302.
[3] Guido Stampacchia (1966), "Équations elliptiques du second ordre à coefficients discontinus" (notes du cours donné à la 4me session du Séminaire de mathématiques supérieures de l'Université de Montréal, tenue l'été 1965), (in French), Séminaire de mathématiques supérieures 16, Montréal: Les Presses de l'Université de Montréal, pp. 326, ISBN 0-8405-0052-1, MR0251373, Zbl 0151.15501.
Best Answer
If $\vec B$ is a constant vector, then we have
$$\begin{align} \nabla \left(\frac{\vec B\cdot \vec r}{r^3}\right)&=\nabla \left(\frac{ B_r}{r^2}\right)\\\\ &=B_r\nabla\left(\frac1{r^2}\right)+\frac1{r^2}\underbrace{\nabla (B_r)}_{=0}\\\\ &=B_r \left(-\frac{2\vec r}{r^4}\right) \end{align}$$
Hence, we have
$$\vec r\times \nabla \left(\frac{\vec B\cdot \vec r}{r^3}\right)=0$$