Reduced row echelon form

Question

Is there a function in R that produces the reduced row echelon form of a matrix?. This reference says there isn't. Do you agree?

Answer 1

I don't have enough rep to comment, but the function given above by soldier.moth in the accepted answer [edit 2018: no longer the accepted answer] is buggy - it doesn't handle matrices where the RREF solution has zeroes on its main diagonal. Try e.g.

m<-matrix(c(1,0,1,0,0,2),byrow=TRUE,nrow=2) rref(m)

and note that the output is not in RREF.

I think I have it working, but you may want to check outputs for yourself:

rref <- function(A, tol=sqrt(.Machine$double.eps),verbose=FALSE,
                 fractions=FALSE){
  ## A: coefficient matrix
  ## tol: tolerance for checking for 0 pivot
  ## verbose: if TRUE, print intermediate steps
  ## fractions: try to express nonintegers as rational numbers
  ## Written by John Fox
  # Modified by Geoffrey Brent 2014-12-17 to fix a bug
  if (fractions) {
    mass <- require(MASS)
    if (!mass) stop("fractions=TRUE needs MASS package")
  }
  if ((!is.matrix(A)) || (!is.numeric(A)))
    stop("argument must be a numeric matrix")
  n <- nrow(A)
  m <- ncol(A)
  x.position<-1
  y.position<-1
  # change loop:
  while((x.position<=m) & (y.position<=n)){
    col <- A[,x.position]
    col[1:n < y.position] <- 0
    # find maximum pivot in current column at or below current row
    which <- which.max(abs(col))
    pivot <- col[which]
    if (abs(pivot) <= tol) x.position<-x.position+1     # check for 0 pivot
    else{
      if (which > y.position) { A[c(y.position,which),]<-A[c(which,y.position),] } # exchange rows
      A[y.position,]<-A[y.position,]/pivot # pivot
      row <-A[y.position,]
      A <- A - outer(A[,x.position],row) # sweep
      A[y.position,]<-row # restore current row
      if (verbose)
        if (fractions) print(fractions(A))
        else print(round(A,round(abs(log(tol,10)))))
      x.position<-x.position+1
      y.position<-y.position+1
    }
  }
  for (i in 1:n)
    if (max(abs(A[i,1:m])) <= tol)
      A[c(i,n),] <- A[c(n,i),] # 0 rows to bottom
  if (fractions) fractions (A)
  else round(A, round(abs(log(tol,10))))
}