Remarks: There are some widely used tricks for this kind of problems. Here is one of them.
Alternative proof:
Using the identity ($q > 0$)
$$\frac{1}{q} = \int_0^\infty \mathrm{e}^{-qt} \,\mathrm{d} t,$$
we have
$$\sum_{i=1}^n\sum_{j=1}^n \frac{1}{x_i + x_j}
= \int_0^\infty \sum_{i=1}^n\sum_{j=1}^n
\mathrm{e}^{-(x_i + x_j)t} \,\mathrm{d} t
= \int_0^\infty (\mathrm{e}^{- tx_1}
+ \cdots + \mathrm{e}^{- tx_n})^2 \,\mathrm{d} t.$$
Similarly, we have
$$\sum_{i=1}^m \sum_{j=1}^m
\frac{1}{y_i + y_j}
= \int_0^\infty (\mathrm{e}^{- ty_1}
+ \cdots + \mathrm{e}^{- ty_m})^2 \,\mathrm{d} t,$$
and
$$\sum_{i=1}^n \sum_{j=1}^m \frac{1}{x_i + y_j}
= \int_0^\infty (\mathrm{e}^{- tx_1}
+ \cdots + \mathrm{e}^{- tx_n})
(\mathrm{e}^{- ty_1}
+ \cdots + \mathrm{e}^{- ty_m}) \,\mathrm{d} t.$$
Thus, we have
$$\mathrm{LHS} - \mathrm{RHS}
= \int_0^\infty (\mathrm{e}^{- tx_1}
+ \cdots + \mathrm{e}^{- tx_n}
- \mathrm{e}^{- ty_1}
- \cdots - \mathrm{e}^{- ty_m})^2 \,\mathrm{d} t \ge 0.$$
We are done.