We now subtract the appropriate multiples of the first row from the rows $i=2,\dots,n$ and the appropriate multiples of the first column from the columns $ji=2,\dots,n$.
Let us first find the multiple of the first row, to be subtracted from the second one, denoting it as ${\rm row}_2$. Since $p_{11}(\lambda) = \lambda$ and $p_{21}(\lambda) = \lambda^2$, we see that ${\rm row}_2$ equals the first row multiplied by $(p_{21}/p_{11})(\lambda) = \lambda$:
$${\rm row}_2 = (p_{21}/p_{11})(\lambda) \begin{bmatrix} \lambda & \lambda^2 \end{bmatrix} = \begin{bmatrix} \lambda^2 & \lambda^3 \end{bmatrix}.$$
We now subtract this from the current second row:
$$\begin{bmatrix} \lambda^2 & \lambda - 1 \end{bmatrix} - {\rm row}_2 = \begin{bmatrix} 0 & -\lambda^3 + \lambda - 1 \end{bmatrix}.$$
Our current $P(\lambda)$ is
$$P(\lambda) = \begin{bmatrix} \lambda & \lambda^2 \\ \color{red}{0} & \color{red}{-\lambda^3 + \lambda - 1} \end{bmatrix}.$$
As before, we subtracting the first column's multiple will now just remove the nonzeroes in the first row.
$$P(\lambda) = \begin{bmatrix} \lambda & \color{red}{0} \\ 0 & \color{red}{-\lambda^3 + \lambda - 1} \end{bmatrix}.$$
Now, this may seem great, but notice that $p_{11} \nmid p_{22}$, so we need further corrections. This is done by adding the second row to the first one, repeating the above procedure. So, we get
$$P(\lambda) = \begin{bmatrix} \color{red}{\lambda} & \color{red}{-\lambda^3 + \lambda - 1} \\ 0 & -\lambda^3 + \lambda - 1 \end{bmatrix} = \begin{bmatrix} \lambda & \color{red}{\lambda \cdot (-\lambda^2 + 1) + (-1)} \\ \color{red}{0} & -\lambda^3 + \lambda - 1 \end{bmatrix}.$$
We see that $q_{12}(\lambda) = -\lambda^2 + 1$ and $r_{12}(\lambda) = 1$. We need to subtract the first column times $q_{12}(\lambda)$ from the second colum:
$$P(\lambda) = \begin{bmatrix} \lambda & \color{red}{-1} \\ 0 & \color{red}{-\lambda^3 + \lambda - 1} \end{bmatrix}.$$
Exchanging these two, we get
$$P(\lambda) = \begin{bmatrix} -1 & \lambda \\ -\lambda^3 + \lambda - 1 & 0 \end{bmatrix}.$$
After eliminating $p_{12}$, we have
$$P(\lambda) = \begin{bmatrix} -1 & \color{red}{0} \\ -\lambda^3 + \lambda - 1 & \color{red}{-\lambda^4 + \lambda^2 - \lambda} \end{bmatrix}.$$
Then we eliminate $p_{21}$:
$$P(\lambda) = \begin{bmatrix} -1 & 0 \\ \color{red}{0} & \color{red}{-\lambda^4 + \lambda^2 - \lambda} \end{bmatrix}.$$
Repeat the process on the bottom right principal submatrix of order $n-1$. In our case, as above, $n = 2$, so the job is almost done, as we have
$$P(\lambda) = \begin{bmatrix} -1 & 0 \\ 0 & -\lambda^4 + \lambda^2 - \lambda \end{bmatrix}.$$
For this to be a Smith form, we need the nonzero diagonal polynomials to be monic (i.e., have a leading coefficient equal to $1$). In our case, we do this by multiplying $P(\lambda)$ by $-1$. In a more generale case, we'd need to multiply it from either left or right by some constant diagonal matrix.
Finally, we obtain
$$P(\lambda) = \begin{bmatrix} 1 & 0 \\ 0 & \lambda^4 - \lambda^2 + \lambda \end{bmatrix}.$$
Best Answer
The computation requires two processes:
Specifically, you are allowed to
The first goal is to reach diagonal form. Let's first work on column 1: using operation 3 for rows repeatedly, you can, by following the Euclidean algorithm, form a row whose first element is the GCD of the elements in column 1. You can then obtain a matrix with the GCD in the $(1,1)$ position and zeroes in the rest of column 1.
Now work on row 1: do the same thing, but using column operations; eventually you will have the GCD of row 1 in the $(1,1)$ position, and zeroes elsewhere in row 1. You will most likely have messed up column 1, but that's OK. Go back and redo column 1, then redo row 1, and repeat until all elements in row and column 1 are 0 except for the $(1,1)$ element. This process is guaranteed to terminate because the GCD gets smaller each time.
Now you can move on to row/column 2, and repeat the process. Continue for row/column 3, and so on, until you have reached diagonal form.
You may not be done at this point, because the diagonal elements may not satisfy the divisibility requirement of the Smith normal form. You can, however, enforce this by some additional moves as in the following example: $$ \begin{bmatrix}8 & 0\\0 & 12\end{bmatrix}\longrightarrow\begin{bmatrix}8 & 8\\0 & 12\end{bmatrix}\longrightarrow\begin{bmatrix}8 & 8\\-8 & 4\end{bmatrix}\longrightarrow\begin{bmatrix}24 & 8\\0 & 4\end{bmatrix}\longrightarrow\begin{bmatrix}24 & 0\\0 & 4\end{bmatrix}\longrightarrow\begin{bmatrix}4 & 0\\0 & 24\end{bmatrix}. $$ The idea is again to use the Eulidean algorithm. After adding column 1 to column 2, you may have to do several row operations to obtain the GCD of the two original diagonal elements. In the example above, only one row operation was needed for this.
Addendum: Details for the particular example in the OP's question. If you add row 2 to row 1, and then multiply row 1 by $-1$, you get a pivot of 1. You can then subtract suitable multiples of row 1 from rows 2 and 4, and end up with $$ \begin{bmatrix} 1 & 561 & -174 & -80 \\ 0 & -3477 & 1080 & 474 \\ 0 & -255 & 81 & 24 \\ 0 & -4182 & 1299 & 570 \end{bmatrix}, $$ which, by subtracting multiples of column 1 from columns 2, 3, and 4, leads to $$ \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & -3477 & 1080 & 474 \\ 0 & -255 & 81 & 24 \\ 0 & -4182 & 1299 & 570 \end{bmatrix}. $$ We next need to do a series of row operations involving rows 2, 3, and 4 that results in the GCD of 3477, 255, and 4182. If we always choose the pivot of smallest nonzero magnitude, the steps would be
You should now have $$ \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 3 & 234 & -1410 \\ 0 & 0 & -651 & 3906 \\ 0 & 0 & -147 & 882 \end{bmatrix}. $$ You can now subtract suitable multiples of column 2 from columns 3 and 4. Then you just have to deal with the $2\times2$ block in the lower right. You should be able to get it from here.