rref.R

rref <- function(A, rows=length(A[,1]), cols=length(A[1,])) {   
  interchange <- function(A, i, j) { row = A[i,]; A[i,] = A[j,]; A[j,] = row; A }
  scale <- function(A, i, c=1) { A[i,] = c*A[i,]; A }
  replace <- function(A, j, i, c=1) { A[j,] = A[j,] + c*A[i,]; A}
  
  for(i in 1:(cols-1)) {
    j = i # only interchange if j >= i
    found = FALSE # look for a non-zero in the col
    while(found==FALSE && j <= rows) { 
      if(A[j,i] != 0) {
        A = scale(A,j,c=1/A[j,i])
        A = interchange(A,i,j)
        found = TRUE
        for(m in 1:rows) { if(m!=i) A = replace(A,m,i,c=-A[m,i]) }
      }
      j = j + 1 
    }
  } 
  A
}

inverse <- function(A, rows=length(A[,1]), cols=length(A[1,])) {
  I = matrix(rep(c(1,rep(0,cols)), rows), nrow=rows, ncol=cols, byrow=TRUE)
  B = matrix(c(A,I), nrow=rows, ncol=cols*2)
  rref(B)[,(cols+1):(2*cols)]
}

multiply <- function(A, B, rows=length(A[,1]), cols=length(B[1,])) {
  AB = matrix(rep(0,rows*cols),nrow=rows,ncol=cols)
  dim = length(A[1,])
  for(i in 1:rows) {
    for(j in 1:cols) {
      for(k in 1:dim) {
        AB[i,j] = AB[i,j] + A[i,k]*B[k,j]
      }
    }
  }
  AB
}

determinant <- function(A, dim=length(A[,1])) {
  if(dim < 2) {
    return(A[0,0])
  } else if(dim == 2) {
    return(A[1,1]*A[2,2] - A[1,2]*A[2,1])
  } else { # recursion
    total = 0
    for(i in 1:dim) {
      total = total + ifelse((i+1)%%2==0, 1, -1)*A[1,i]*det(A[1:dim!=1, 1:dim!=i])
    }
    return(total)
  }
}

eigenvalue <- function(A, rows=length(A[,1]), cols=length(A[1,])) {
  newtons <- function(f, x, n=5, epsilon=0.00001) {
    for(i in 1:n) {
      x = x - (f(x+epsilon)-f(x))/(epsilon*f(x))
    }
    x
  }
  
  lambdaI <- function(lambda=0) {
    det(A - matrix(rep(c(lambda,rep(0,cols)), rows), nrow=rows, ncol=cols, byrow=TRUE))
  }
  newtons(lambdaI, 0)
}

A = matrix(c(c(0,-5,8,3),c(-2,-7,3,4),c(4,2,7,5)),nrow=3, ncol=4,byrow=TRUE)
rref(A) # demo rref

##      [,1] [,2]       [,3]    [,4]
## [1,]    1    0  8.882e-16  3.9355
## [2,]    0    1 -2.220e-16 -2.0968
## [3,]    0    0  1.000e+00 -0.9355

# 
# A = matrix(c(c(0,1,2),c(1,0,3),c(4,-3,8)),nrow=3,ncol=3,byrow=TRUE)
# inverse(A) # demo inverse
# 
# A1 = matrix(c(3,1,5,2,1,4,3,2,3), nrow=3, ncol=3, byrow=TRUE)
# B1 = matrix(c(1,4,1,2,3,1,5,1,4), nrow=3, ncol=3, byrow=TRUE) 
# multiply(A1,B1) # demo multiply
# 
# A = matrix( c(c(3.0, 6.0, 7.0, 6.0, 3.0, 6.0, 3.0, 6.0),  c(5.0, 3.0, 6.0, 3.0, 1.0, 5.0, 3.0, 8.0),
#               c(2.0, 5.0, 7.0, 3.0, 5.0, 2.0, 8.0, 5.0),  c(1.0, 0.0, 7.0, 5.0, 3.0, 0.0, 6.0, 5.0),  
#               c(4.0, 6.0, 4.0, 2.0, 7.0, 8.0, 4.0, 3.0),  c(3.0, 9.0, 5.0, 3.0, 2.0, 1.0, 5.0, 7.0),  
#               c(3.0, 8.0, 5.0, 3.0, 8.0, 0.0, 9.0, 7.0),  c(5.0, 2.0, 1.0, 2.0, 3.0, 5.0, 5.0, 3.0)),
#             nrow=8, ncol=8, byrow=TRUE)
# det(A)
# A = matrix(c(c(5,-2,6,-1),c(0,3,-8,0),c(0,0,5,4),c(0,0,0,1)), nrow=4, ncol=4, byrow=TRUE)
# eigen(A)