Reduced Row-Echelon Form of a Matrix

Question

The goal of this challenge is to create a program that takes in a matrix and outputs its reduced row-echelon form.

A matrix is in reduced row-echelon form if it meets all of the following conditions:

1. If there is a row where every entry is zero, then this row lies below any other row that contains a nonzero entry.
2. The leftmost nonzero entry of a row is equal to 1.
3. The leftmost nonzero entry of a row is the only nonzero entry in its column.
4. Consider any two different leftmost nonzero entries, one located in row i, column j and the other located in row s, column t. If s>i, then t>j.

Source

The general process to transform the matrix is:

1. Deal with each row i from 1 to n in turn, and work across the columns j from 1 to m skipping any column of all zero entries.
2. Find the next column j with a nonzero entry.
3. Interchange rows, if necessary, so that the pivot element A(i,j) is nonzero.
4. Make the pivot equal to 1 by dividing each element in the pivot row by the value of the pivot.
5. Make all elements above and below the pivot equal to 0 by subtracting a suitable multiple of the pivot row from each other row.
6. Repeat for all rows.

If you want to read more on this type of matrix, view the Wikipedia article on it and a tool and article(steps above) that shows the steps to transform the matrix.

As for the actual challenge, here it goes:

The input can be given in any way you wish through STDIN or equivalent, just please explain it in your answer. The output will be the reduced row-echelon form of the input in the same form as the input through STDOUT or equivalent. Standard loopholes are not allowed and external libraries or functions that do this task are also not allowed (TI-BASIC's rref( command, for example). You may write a full program or function. This is code golf, lowest BYTES wins. Good luck!

Example Input: [[2,1,1,14][-1,-3,2,-2][4,-6,3,-5]]

Example Output: [[1,0,0,1][0,1,0,5][0,0,1,7]]

Answer 1

R, 232 bytes

function(M){p=1;R=nrow(M);C=ncol(M);for(r in 1:R){if(C<=p)break;i=r;while(M[i,p]==0){i=i+1;if(R==i){i=r;p=p+1;if(C==p)return(M)}};t=M[i,];M[i,]=M[r,];M[r,]=t;M[r,]=M[r,]/M[r,p];for(i in 1:R)if(i!=r)M[i,]=M[i,]-M[r,]*M[i,p];p=p+1};M}

This is just an implementation of the usual Gaussian elimination algorithm.

Ungolfed:

rref <- function(M) {
    p <- 1
    R <- nrow(M)
    C <- ncol(M)
    for (r in 1:R) {
        if (C <= p)
            break
        i <- r
        while (M[i, p] == 0) {
            i <- i + 1
            if (R == i) {
                i <- r
                p <- p + 1
                if (C == p)
                    return(M)
            }
        }
        t <- M[i, ]
        M[i, ] <- M[r, ]
        M[r, ] <- t
        M[r, ] <- M[r, ] / M[r, p]
        for (i in 1:R)
            if (i != r)
                M[i, ] <- M[i, ] - M[r, ] * M[i, p]
        p <- p + 1
    }
    M
}