Generate random mumbers with bivariate gamma distribuition

How to generate random numbers with bivariate gamma distribuition. The density is:

$$f_{(X,Y)}(x,y) = \alpha^{p+q} x^{p-1} (y-x)^{q-1} \exp{-\alpha y} / [\Gamma(p) \Gamma(q)]\,\mathbb{I}_{0\le x\le y}$$

with $y>x>0,\ \alpha>0,\ p>0$ and $q>0$.

I did not find any package on R that do this and nothing in literature.