I want to calculate the matrix elements of the operator $\hat{x} \hat{p}$ in momentum and position basis, that is the two quantities ($p,q$ - momenta, $x,y$ - positions):
$$\langle p|\hat{x} \hat{p}|q\rangle$$ $$\langle x|\hat{x} \hat{p}|y\rangle$$
I don't know how to do this. I write $\hat{p}|q\rangle = q | q \rangle$. And $\hat{x} |q \rangle = -i\hbar \frac{d}{dp} | q\rangle$, so
$$\langle p|\hat{x} \hat{p}|q\rangle = -i\hbar q \frac{d}{dp} \delta(p-q)$$
This is nonsensical.
How do I proceed?
It's not non-sensical at all, except that there shouldn't be a minus sign (as mentioned in the comments) and that you took an operator outside of an expectation value, which I think worked out OK in this case but in general should be avoided. More conservatively,
$$ \hat x = i \hbar \frac{d}{d \hat p} $$
it follows that
$$ \langle q \mid \hat x \hat p \mid q' \rangle = q' \langle q \mid \hat x \mid q' \rangle = i \hbar q' \langle q \mid \frac{d}{d \hat p} \mid q' \rangle = i \hbar q' \langle q \mid \frac{d}{d q'} \mid q' \rangle $$
Remember, an actual physical state is never an idealized momentum (or position) eigenket. A slightly more realistic description of a physical state is a wave packet with a narrow width around some momentum:
$$ \mid p, \sigma \rangle = \int \frac{dq}{2 \pi \hbar} \, \left( \frac{2 \pi \hbar^2}{\sigma^2}\right)^{1/4} e^{- \frac{1}{4} \left( \frac{q-p}{\sigma} \right)^2} \mid q \rangle $$
You can check that this is normalized to unity, because
$$ \langle p, \sigma \mid p, \sigma \rangle = 1 $$
Remember that $\langle p \mid q \rangle = 2 \pi \hbar \delta(p-q)$. Now, the expectation value of $\hat x \hat p$ for this physical state is
$$ \langle p, \sigma \mid \hat x \hat p \mid p, \sigma \rangle = \frac{1}{\sqrt{2 \pi} \, \sigma } \int dq \, dq' \, e^{- \frac{1}{4} \left( \frac{q-p}{\sigma} \right)^2} e^{- \frac{1}{4} \left( \frac{q'-p}{\sigma} \right)^2} \frac{\langle q \mid \hat x \hat p \mid q' \rangle}{2 \pi \hbar} $$
Using the above result, we can integrate the $q'$ integral by parts, yielding
$$ \frac{-i\hbar}{\sqrt{2 \pi} \, \sigma } \int dq \, dq' \, e^{- \frac{1}{4} \left( \frac{q-p}{\sigma} \right)^2} \frac{d}{dq'} \left( q' e^{- \frac{1}{4} \left( \frac{q'-p}{\sigma} \right)^2} \right) \frac{\langle q \mid q' \rangle}{2 \pi \hbar} $$
Since $\langle q \mid q' \rangle = 2 \pi \hbar \delta(q-q')$, we can integrate out q (and then relabel q' back to q), making this a single integral
$$ \frac{-i\hbar}{\sqrt{2 \pi} \, \sigma} \int dq \, e^{- \frac{1}{4} \left( \frac{q-p}{\sigma} \right)^2} \frac{d}{dq} \left( q e^{- \frac{1}{4} \left( \frac{q-p}{\sigma} \right)^2} \right) $$
We can integrate by parts again, ending up with
$$ \frac{i\hbar}{\sqrt{2 \pi} \, \sigma} \int dq \, q e^{- \frac{1}{4} \left( \frac{q-p}{\sigma} \right)^2} \frac{d}{dq} \left( e^{- \frac{1}{4} \left( \frac{q-p}{\sigma} \right)^2} \right) $$
which is a bit easier to evaluate since we don't have to use the product rule. The derivative brings down a factor $-(q-p)/2\sigma^2$, and the two exponentials combine, leaving us with
$$ \frac{-i\hbar}{ 2\sqrt{2 \pi} \, \sigma^3} \int dq \, q(q-p) e^{- \frac{1}{2} \left( \frac{q-p}{\sigma} \right)^2} $$
The first $q$ in the integrand can be rewritten as $(q-p) + p$. Then there are two terms, the latter of which has an odd integrand (of the form $(q-p) \exp(\alpha (q-p)^2)$ and vanishes. So, changing the integrand to $x = (q-p)/\sigma$, we are left with
$$ \frac{-i\hbar}{ 2\sqrt{2 \pi}} \int dx \, x^2 e^{- \frac{1}{2} x^2} $$
The remaining integral is a textbook Gaussian integral, equal to
$$ \int dx \, x^2 e^{-\frac{1}{2} x^2} = \sqrt{2\pi} $$
So, the expectation value of the operator $\hat x \hat p$ for the wave packet is simply
$$ \frac{-i\hbar}{2} $$
