Puedes seguir el algoritmo que da la Wikipedia
aquí, y sale.
En efecto la idea es usar la identidad de Bézout. Tenemos que (uso la misma notación que en la Wikipedia) \( \sigma a + \tau b = d \), de donde \( \sigma \alpha + \tau \gamma = 1 \), con \( \alpha := a/d \) y \( \gamma := b/d \). Entonces la matriz \( \begin{pmatrix}{\sigma}&{\tau}\\{-\gamma}&{\alpha}\end{pmatrix} \) es invertible (porque tiene determinante \( 1 \)), y si multiplicas por la izquierda con la del problema te queda:
\[ \begin{pmatrix}{d}&{-4\sigma b + \tau a}\\{0}&{\frac{a^2+4b^2}{d}}\end{pmatrix} \].
Y ahora ya estás, porque la entrada de fuera de la diagonal es un múltiplo de \( d \) (al serlo tanto \( a \) como \( b \)), y por tanto restando de la segunda columna un múltiplo adecuado de la primera columna (que corresponde a multiplicar por una matriz invertible por la derecha) puedes hacer cero la entrada de fuera de la diagonal. Y ya te queda la matriz que te dan.