Your manual solution, if carried further, yields
Solve[r[th] + Log[r[th]] == th + Log[th] + E, r[th]]
(* {{r[th] -> ProductLog[E^(E + th) th]}} *)
which is equivalent to the DSolve
solution.
You can also verify the DSolve
solution with
{y'[x] == (x y[x] + y[x])/(x y[x] + x), y[1] == E} /. msol // FullSimplify // N[#, 80] &
(For some reason Mathematica does not simplify ProductLog[E^(1 + E)]
, so I verified it numerically out to 80 digits -- an arbitrary choice, but hopefully convincing.)