Here is a Scheme program to compute the Reduced Row Echelon Form of a given matrix. Scheme can do exact rational arithmetic and this circumvents the difficult issue of threshold tests for being ‘close enough to zero’ to be zero. With exact input the output is also exact.

See a geometric interpretation of this idea. Here is a repository version.

The reductions allowed here are:

If each row of a matrix is taken as the coordinates in some basis of some vector space V, then these reductions leave unchanged the sub-space of V spanned by the set of such vectors. Any matrix can be reduced by a sequence of such reductions, to a matrix in RREF and that matrix is unique.

The argument, x, and the yield of (rref x), are each a list of rows, each of which is a list of numbers. (rank x) returns the rank of a matrix x. Viewing x as a list of vectors, (minspan x) returns a minimum sublist of x spanning what x spans. Members of x that are linear combinations of earlier members are omitted from the yield.

x and (rref x), viewed each as a list of vectors relative to some constant basis, span the same subspace. This is the ‘alias’ perspective. (length x) = (length (rref x)). Removing trailing zero vectors (rows) from a matrix in RREF results in a canonical spanning set for the space that the remaining vectors span. This is canonical only while the basis remains constant. Arbitrary subspaces can be defined by vector spanning sets. Such a clipped matrix might be used to denote a Grassmann manifold value. 

(define (rref x) (if (or (null? x) (null? (car x))) x
  (let ((a (call-with-current-continuation (lambda (cc) 
    (let df ((y x)) (if (null? y) (cc x) (if (zero? (caar y)) (let ((z (df (cdr y))))
       (if (null? z) y (cons (car z) (cons (car y)(cdr z))))) y)))))))
  (if (zero? (caar a)) (map (lambda (x) (cons 0 x)) (rref (map cdr a))) (let* (
      (b (let* ((r (/ (caar a))) (top (map (lambda (z) (* z r)) (car a))))
        (cons top (map (lambda (row) (map (lambda (t re) (- re (* t (car row))))
          top row)) (cdr a)))))
      (c (rref (map cdr (cdr b))))
  ; now we must ensure a zero in row (cdar b) above each initial 1 a row of c.
      (d (let w ((sm (cdar b))(i c))
        (if (null? i) sm (w (let sr ((x sm) (y (car i)))
          (if (null? x) '() (if (zero? (car y)) (cons (car x) (sr (cdr x)(cdr y)))
        (map (lambda (p q)(- p (* (car x) q))) x y)))) (cdr i))))))
   (cons (cons 1 d) (map (lambda (x) (cons 0 x)) c)))))))

(define (rank x)(let ex ((w (rref x))(n 0))
  (if (or (null? w) (let ze ((p (car w))) (or (null? p)
      (and (zero? (car p)) (ze (cdr p)))))) n
   (ex (cdr w)(+ n 1)))))

(define (minspan x) (let nx ((a (rref (transpose x)))) (if (null? a) '()
  (let ((b (let cr ((x (car a))(y x)) (if (null? x) '()
  (if (zero? (car x)) (cr (cdr x)(cdr y)) (car y))))))
    (if (null? b) '() (cons b (nx (cdr a))))))))
The following five functions are the form of the first seemingly correct version of rref. They may thus help explain the code above. 
; prm takes a matrix and moves (if it can) to the top, a row that doesn't start with 0.
(define (prm a) (if (or (null? a) (null? (car a))) '()
   (let ((a (call-with-current-continuation (lambda (cc) 
  (let df ((y a)) (if (null? y) (cc a) (if (not (zero? (caar y))) y (let ((z (df (cdr y))))
    (if (null? z) y (cons (car z) (cons (car y)(cdr z)))))))))))) a)))
(prm '((0 2)(1 4)(0 5))) ;=> ((1 4) (0 2) (0 5))
; =======
; zz takes a matrix x where (caar x) isn’t 0 and makes left column (1 0 0 ...).
(define (zz x) (let* ((r (/ (caar x))) (top (map (lambda (z) (* z r)) (car x))))
  (cons top (map (lambda (row) (map (lambda (t re) (- re (* t (car row))))
     top row)) (cdr x)))))
; =======
; yy takes a matrix with a left column of 0's and performs rref on all but left row.
(define (yy a) (map (lambda (x) (cons 0 x)) (rref (map cdr a))))
; =======
; sr takes rows x and y. Initial non zero element of y, if any, is 1.
; sr returns row z which is x+py such that element in z above 1 in y is 0.
(define (sr x y) (if (null? x) '() (if (zero? (car y)) (cons (car x) (sr (cdr x)(cdr y)))
   (map (lambda (p q)(- p (* (car x) q))) x y))))
; (sr '(3 4 5 6 7 8) '(0 0 1 2 3 4)) ; => (3 4 0 -4 -8 -12)
; =======
; rref takes a matrix x and returns its RREF form.
(define (rref x) (if (or (null? x) (null? (car x))) x
  (let ((a (prm x))) (if (zero? (caar a)) (yy a)
  (let* ((b (zz a)) (c (rref (map cdr (cdr b))))
  ; now we must ensure a zero in row (cdar b) above each initial 1 a row of c.
    (d (let w ((sm (cdar b))(i c)) (if (null? i) sm (w (sr sm (car i)) (cdr i))))))
   (cons (cons 1 d) (map (lambda (x) (cons 0 x)) c)))))))
The C code below generates a random matrix with known RREF. 
#include <stdio.h>
#include < stdlib.h>
int main(){
int z[][10] = {{0, 1, 2, 6, 0, 1, 0, 0, 4, 0}, {0, 0, 0, 0, 1, 1, 0, 0, 1, 1},
 {0, 0, 0, 0, 0, 0, 1, 0, 4, 1}, {0, 0, 0, 0, 0, 0, 0, 1, 2, 1},
 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};

void pr() {int j; printf("("/*)*/); for(j=0; j<5; ++j) {int k; printf("("/*)*/);
  for(k=0; k<10; ++k) printf("%d ", z[j][k]); printf(/*(*/")\n");}
    printf(/*(*/")\n");}
pr();
{int k = 30; while(k--) {
  int i = random()%5, j= random()%5, n=10; int s = 1+random()%3;
  while(n--) z[i][n] += s*z[j][n];}}
pr(); return 0;}

yields
((0 1 2 6 0 1 0 0 4 0 )
 (0 0 0 0 1 1 0 0 1 1 )
 (0 0 0 0 0 0 1 0 4 1 )
 (0 0 0 0 0 0 0 1 2 1 )
 (0 0 0 0 0 0 0 0 0 0 )
)
((0 1176 2352 7056 3032 4208 2160 2376 21128 7568 )
(0 169 338 1014 435 604 310 341 3033 1086 )
(0 4068 8136 24408 10494 14562 7488 8220 73158 26202 )
(0 2852 5704 17112 7355 10207 5246 5762 51271 18363 )
(0 6 12 36 15 21 10 12 103 37 )
)
Another check: (rref '((3 4 5 6 7 8)(0 0 1 2 3 4))) ;=> ((1 4/3 0 -4/3 -8/3 -4) (0 0 1 2 3 4))