Gröbner basis of $f=x^2 y +xy^2 +y^2$, $f_1 =y^2 -1$

Consider the following polynomials:

$$f=x^2 y +xy^2 +y^2; \qquad f_1 =y^2 -1; \quad\text{and}\quad f_2 = xy -1.$$

Dividing by $(f_1,f_2)$, one gets

$$f=f_1 (x+1)+f_2 (x)+2x+1.$$

During lecture the teacher found a Gröbner basis by multiplying 3 matrices:

$$\begin{bmatrix} -1 &-1& -y \\ 0& 1& 0 \\ 0& 0& 1\end{bmatrix} \begin{bmatrix} -1 &0 \\ 0& 1\\ y& -x\end{bmatrix} \begin{bmatrix} xy -1\\ y^2-1\end{bmatrix} = \begin{bmatrix} 0 \\ y^2-1\\ x-y \end{bmatrix}$$

This is supposed to be a Gröbner basis but I'm not sure how it was obtained through matrix multiplication.

This is an example from Cox, Little, and O'Shea's, Ideals, Varieties, and Algorithms, chapter 2, section 3, example 4.

We have here a case where the Buchberger algorithm is cast into a linear algebra system. It is especially applied to a system of two polynomials. So, what do these matrices represent? The second matrix transforms a set of two equations to a set of three; one with changed sign: $-f_2$, one identical $f_1$ and a new one as combination $\alpha f_1 + \beta f_2$. Whence the following setup:
$$\begin{pmatrix} -1 & 0\\ 0 & 1 \\ \beta & \alpha \end{pmatrix} \begin{pmatrix} f_2 \\f_1\end{pmatrix} = \begin{pmatrix} -f_2\\ f_1 \\ \alpha f_1 + \beta f_2 \end{pmatrix}.$$
The coefficients $\alpha$ and $\beta$ are determined by the Buchbeger algorithm and construct a so-called $S$-polynomial. More details here. The $S$-polynomial is defined by (as done in the link): $$\operatorname{SPOL}(f_1,f_2) = \frac{\operatorname{LCM}(\operatorname{LT}(f_1),\operatorname{LT}(f_2))}{\operatorname{LM}(f_1)}f_1-\frac{\operatorname{LCM}(\operatorname{LT}(f_1),\operatorname{LT}(f_2))}{\operatorname{LM}(f_2)}f_2,$$ where $\operatorname{LCM}$ means the least common multiple, $\operatorname{LT}(f)$ the leading term of $f$ and $\operatorname{LM}(f)$ the leading monomial of $f$. This gives us the following values for $\alpha$ and $\beta$: $$\alpha = \frac{\operatorname{LCM}(\operatorname{LT}(f_1),\operatorname{LT}(f_2))}{\operatorname{LM}(f_1)},\ \beta = -\frac{\operatorname{LCM}(\operatorname{LT}(f_1),\operatorname{LT}(f_2))}{\operatorname{LM}(f_2)}.$$ So in the example of the question $\alpha = \frac{xy^2}{y^2} = x$ and $\beta = -\frac{xy^2}{xy}=-y$. Now $$\operatorname{SPOL}(f_1,f_2) = \alpha f_1 + \beta f_2 = x(y^2-1) - y(xy-1) = y - x = -\operatorname{SPOL}(f_2,f_1),$$ so the product of the last two matrices gives us: $$\begin{pmatrix} -f_1 \\ f_2 \\ \operatorname{SPOL}(f_2,f_1) \end{pmatrix} .$$ The significance of the first matrix is the decision which of the two original polynomials remains and which will be replaced by $0$. Let us call our new element of the Gröbner basis $f_3$ and let's calculate $\operatorname{SPOL}(f_2,f_3) = \frac{xy}{xy}(xy-1)-\frac{xy}{x}(x-y)=y^2-1 = f_1$. Hence the first row in the first matrix. This proves that if we omit $f_2$ the remaining polynomials form a Gröbner basis. If on the other hand we calculate $\operatorname{SPOL}(f_1,f_3) = -f_3$ so nothing new. The result is that we have obtained a Gröbner basis simply by multiplying matrices.