In [110]: 
from sympy import *
In [111]: 
A=Matrix([[1,0,1,0,1,0,3],
          [0,1,0,1,0,1,2],
          [1,1,1,-1,-1,-1,1],
          [1,0,1,0,-1,-1,4],
          [0,0,0,0,1,2,0],
          [0,0,0,0,1,-1,-1]])
A
Out[111]:
Matrix([
[1, 0, 1,  0,  1,  0,  3],
[0, 1, 0,  1,  0,  1,  2],
[1, 1, 1, -1, -1, -1,  1],
[1, 0, 1,  0, -1, -1,  4],
[0, 0, 0,  0,  1,  2,  0],
[0, 0, 0,  0,  1, -1, -1]])
In [114]: 
B=A.rref()[0];
B #B is the reduced row echelon form of A
Out[114]:
Matrix([
[1, 0, 1, 0, 0, 0, 11/3],
[0, 1, 0, 0, 0, 0, -2/3],
[0, 0, 0, 1, 0, 0,  7/3],
[0, 0, 0, 0, 1, 0, -2/3],
[0, 0, 0, 0, 0, 1,  1/3],
[0, 0, 0, 0, 0, 0,    0]])
In [119]: 
#Next we're going to reproduce B using scales and skews
#It will take about 20 steps.
A1=A
A1
Out[119]:
Matrix([
[1, 0, 1,  0,  1,  0,  3],
[0, 1, 0,  1,  0,  1,  2],
[1, 1, 1, -1, -1, -1,  1],
[1, 0, 1,  0, -1, -1,  4],
[0, 0, 0,  0,  1,  2,  0],
[0, 0, 0,  0,  1, -1, -1]])
In [120]: 
def row_op(M,a,i,b=0,j=1): 
    '''Returns a new matrix with
    row i of M replaced by a*(row i)+b*(row j),
    or by just a*(row i) if b, j omitted.'''
    N=Matrix(M)
    N.zip_row_op(i-1,j-1,lambda x,y: a*x+b*y)
    return N
In [121]: 
A2=row_op(A1,1,3,-1,1) #R_3 --> R_3 - R_1
A3=row_op(A2,1,4,-1,1) #R_4 --> R_4 - R_1
A3
Out[121]:
Matrix([
[1, 0, 1,  0,  1,  0,  3],
[0, 1, 0,  1,  0,  1,  2],
[0, 1, 0, -1, -2, -1, -2],
[0, 0, 0,  0, -2, -1,  1],
[0, 0, 0,  0,  1,  2,  0],
[0, 0, 0,  0,  1, -1, -1]])
In [122]: 
A4=row_op(A3,1,3,-1,2) #R_3 --> R_3 - R_2
A4
Out[122]:
Matrix([
[1, 0, 1,  0,  1,  0,  3],
[0, 1, 0,  1,  0,  1,  2],
[0, 0, 0, -2, -2, -2, -4],
[0, 0, 0,  0, -2, -1,  1],
[0, 0, 0,  0,  1,  2,  0],
[0, 0, 0,  0,  1, -1, -1]])
In [123]: 
A5=row_op(A4,2,2,1,3) #R_2 --> 2R_2 + R_3
A5
Out[123]:
Matrix([
[1, 0, 1,  0,  1,  0,  3],
[0, 2, 0,  0, -2,  0,  0],
[0, 0, 0, -2, -2, -2, -4],
[0, 0, 0,  0, -2, -1,  1],
[0, 0, 0,  0,  1,  2,  0],
[0, 0, 0,  0,  1, -1, -1]])
In [124]: 
A6=row_op(A5,2,1,1,4)  #R_1 --> 2R_1 + R_4
A7=row_op(A6,1,2,-1,4) #R_2 -->  R_2 - R_4
A8=row_op(A7,1,3,-1,4) #R_3 -->  R_3 - R_4
A9=row_op(A8,2,5,1,4)  #R_5 --> 2R_5 + R_4
A10=row_op(A9,2,6,1,4) #R_6 --> 2R_6 + R_4
A10
Out[124]:
Matrix([
[2, 0, 2,  0,  0, -1,  7],
[0, 2, 0,  0,  0,  1, -1],
[0, 0, 0, -2,  0, -1, -5],
[0, 0, 0,  0, -2, -1,  1],
[0, 0, 0,  0,  0,  3,  1],
[0, 0, 0,  0,  0, -3, -1]])
In [125]: 
A11=row_op(A10,3,1, 1,5) #R_1 --> 3R_1 + R_5
A12=row_op(A11,3,2,-1,5) #R_2 --> 3R_2 - R_5
A13=row_op(A12,3,3, 1,5) #R_3 --> 3R_3 + R_5
A14=row_op(A13,3,4, 1,5) #R_4 --> 3R_4 + R_5
A15=row_op(A14,1,6, 1,5) #R_6 -->  R_6 + R_5
A15
Out[125]:
Matrix([
[6, 0, 6,  0,  0, 0,  22],
[0, 6, 0,  0,  0, 0,  -4],
[0, 0, 0, -6,  0, 0, -14],
[0, 0, 0,  0, -6, 0,   4],
[0, 0, 0,  0,  0, 3,   1],
[0, 0, 0,  0,  0, 0,   0]])
In [127]: 
A16=row_op(A15, Rational(1,6),1) #R_1 -->  R_1/6
A17=row_op(A16, Rational(1,6),2) #R_2 -->  R_2/6
A18=row_op(A17,-Rational(1,6),3) #R_3 --> -R_3/6
A19=row_op(A18,-Rational(1,6),4) #R_4 --> -R_4/6
A20=row_op(A19, Rational(1,3),5) #R_5 -->  R_4/3
A20
Out[127]:
Matrix([
[1, 0, 1, 0, 0, 0, 11/3],
[0, 1, 0, 0, 0, 0, -2/3],
[0, 0, 0, 1, 0, 0,  7/3],
[0, 0, 0, 0, 1, 0, -2/3],
[0, 0, 0, 0, 0, 1,  1/3],
[0, 0, 0, 0, 0, 0,    0]])
In [128]: 
B #exactly the same as A20
Out[128]:
Matrix([
[1, 0, 1, 0, 0, 0, 11/3],
[0, 1, 0, 0, 0, 0, -2/3],
[0, 0, 0, 1, 0, 0,  7/3],
[0, 0, 0, 0, 1, 0, -2/3],
[0, 0, 0, 0, 0, 1,  1/3],
[0, 0, 0, 0, 0, 0,    0]])
In []: