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 []: