Incomplete Gauss Jordan elimination: what have I left out

Question

I wrote some C++ to implement a Gauss-Jordan elimination: it's all a bit rough and ready (I wanted to get a fix on the algorithm rather than write a piece of "production" code) but you can see it here: https://github.com/mcmenaminadrian/nocsim/blob/master/gjsolver.cpp

When handling this set of data:

variables: -1,2,-2,0,-2,-2,0,0,0,2

coefficients (the last in each line is the sum of the equation, 
i.e., it's a number not a coefficient as such): 

-2,0,-2,-1,1,-2,1,-2,1,-1,6
1,0,0,2,1,0,1,2,0,-2,-7
0,-2,0,2,1,-2,1,2,0,1,0
0,0,2,-2,1,2,-1,-2,1,1,-8
1,2,0,-2,1,-1,1,-2,-1,-1,1
0,-1,-1,-1,-2,0,0,0,0,2,8
2,1,1,1,-2,2,-2,2,0,2,2
0,-1,-1,1,1,0,1,1,1,1,0
-2,-2,1,-1,-1,1,1,0,0,0,-4
-1,-1,-1,-1,2,-1,-1,-1,0,0,-1

The solver outputs this:

DIAGONAL FORM
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,  == -1 / 1
0 , -967/283 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,  == -1934 / 283
0 , 29/92 , 283/92 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,  == -127 / 23
0 , -79/46 , 19/46 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,  == -98 / 23
0 , -9/23 , 1/23 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,  == -66 / 23
0 , -25/46 , 59/46 , 0 , 0 , 1 , 0 , 0 , 0 , 0 ,  == -130 / 23
0 , -9/23 , 1/23 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ,  == -20 / 23
0 , 59/46 , -27/46 , 0 , 0 , 0 , 0 , 1 , 0 , 0 ,  == 86 / 23
0 , 181/46 , -61/46 , 0 , 0 , 0 , 0 , 0 , 2 , 0 ,  == 242 / 23
0 , 5/23 , -21/23 , 0 , 0 , 0 , 0 , 0 , 1 , 1 ,  == 98 / 23

Which looks correct to me - but I can also see, by simple visual inspection, that we can solve this fully for all variables, as once we work out the second variable fully, the rest fall into place.

My question (out of interest rather than any pressing need), is what steps am I missing in the simple algorithm (set out below) that I am using:

1. Attempt to create upper triangular matrix in standard manner, but with no row re-ordering (so if coefficient in pivot is 0, then move on to the next row/diagonal point)

2. Attempt to eliminate non-diagonal elements column by column in upper triangle by subtracting appropriate multiples of pivot from rows above (and where, as a result of 0 diagonals in step 1 above, also subtracting the multiple of the residual elements in the lower triangle).

Now these two steps work perfectly for a matrix where none of the diagonal/pivot elements are 0, and also seem to give a correct, if not full, result for a matrix where an element is 0, but it looks as though I am missing a further step - is it just a question of renormalising the row so the diagonal element is 1 before attempting step 2?