Here is an example of Gauss-Jordan Elimination with pari. Each elimination step is accomplished by a row operation which is interpreted as matrix multiplication on the left by an elementary matrix. For example, the first operation is replacing the second row by itself plus -2 times the first. Check that E1 has the correct effect. (In this text version the columns in the output matrices don't line up very nicely.) (16:01) gp > A=[1,3,-2,0,2,0,0;2,6,-5,-2,4,-3,-1;0,0,5,10,0,15,5;2,6,0,8,4,18,6] %37 = [1 3 -2 0 2 0 0] [2 6 -5 -2 4 -3 -1] [0 0 5 10 0 15 5] [2 6 0 8 4 18 6] (16:03) gp > E1=[1,0,0,0;-2,1,0,0;0,0,1,0;0,0,0,1] %38 = [1 0 0 0] [-2 1 0 0] [0 0 1 0] [0 0 0 1] (16:04) gp > A1=E1*A %39 = [1 3 -2 0 2 0 0] [0 0 -1 -2 0 -3 -1] [0 0 5 10 0 15 5] [2 6 0 8 4 18 6] (16:04) gp > E2=[1,0,0,0;0,1,0,0;0,0,1,0;-2,0,0,1] %40 = [1 0 0 0] [0 1 0 0] [0 0 1 0] [-2 0 0 1] (16:05) gp > A2=E2*A1 %41 = [1 3 -2 0 2 0 0] [0 0 -1 -2 0 -3 -1] [0 0 5 10 0 15 5] [0 0 4 8 0 18 6] (16:05) gp > E3=[1,0,0,0;0,-1,0,0;0,0,1,0;0,0,0,1] %42 = [1 0 0 0] [0 -1 0 0] [0 0 1 0] [0 0 0 1] (16:06) gp > A3=E3*A2 %43 = [1 3 -2 0 2 0 0] [0 0 1 2 0 3 1] [0 0 5 10 0 15 5] [0 0 4 8 0 18 6] (16:06) gp > E4=[1,2,0,0;0,1,0,0;0,0,1,0;0,0,0,1] %44 = [1 2 0 0] [0 1 0 0] [0 0 1 0] [0 0 0 1] (16:06) gp > A4=E4*A3 %45 = [1 3 0 4 2 6 2] [0 0 1 2 0 3 1] [0 0 5 10 0 15 5] [0 0 4 8 0 18 6] (16:06) gp > E5=[1,0,0,0;0,1,0,0;0,-5,1,0;0,0,0,1] %46 = [1 0 0 0] [0 1 0 0] [0 -5 1 0] [0 0 0 1] (16:07) gp > A5=E5*A4 %47 = [1 3 0 4 2 6 2] [0 0 1 2 0 3 1] [0 0 0 0 0 0 0] [0 0 4 8 0 18 6] (16:07) gp > E6=[1,0,0,0;0,1,0,0;0,0,1,0;0,-4,0,1] %48 = [1 0 0 0] [0 1 0 0] [0 0 1 0] [0 -4 0 1] (16:08) gp > A6=E6*A5 %49 = [1 3 0 4 2 6 2] [0 0 1 2 0 3 1] [0 0 0 0 0 0 0] [0 0 0 0 0 6 2] (16:08) gp > E7=[1,0,0,0;0,1,0,0;0,0,0,1;0,0,1,0] %50 = [1 0 0 0] [0 1 0 0] [0 0 0 1] [0 0 1 0] (16:08) gp > A7=E7*A6 %51 = [1 3 0 4 2 6 2] [0 0 1 2 0 3 1] [0 0 0 0 0 6 2] [0 0 0 0 0 0 0] (16:08) gp > E8=[1,0,0,0;0,1,0,0;0,0,1/6,0;0,0,0,1] %52 = [1 0 0 0] [0 1 0 0] [0 0 1/6 0] [0 0 0 1] (16:09) gp > A8=E8*A7 %53 = [1 3 0 4 2 6 2] [0 0 1 2 0 3 1] [0 0 0 0 0 1 1/3] [0 0 0 0 0 0 0] (16:09) gp > E9=[1,0,0,0;0,1,-3,0;0,0,1,0;0,0,0,1] %54 = [1 0 0 0] [0 1 -3 0] [0 0 1 0] [0 0 0 1] (16:10) gp > A9=E9*A8 %55 = [1 3 0 4 2 6 2] [0 0 1 2 0 0 0] [0 0 0 0 0 1 1/3] [0 0 0 0 0 0 0] (16:11) gp > E10=[1,0,-6,0;0,1,0,0;0,0,1,0;0,0,0,1] %58 = [1 0 -6 0] [0 1 0 0] [0 0 1 0] [0 0 0 1] (16:11) gp > A10=E10*A9 %59 = [1 3 0 4 2 0 0] [0 0 1 2 0 0 0] [0 0 0 0 0 1 1/3] [0 0 0 0 0 0 0]