Matrix LU Factorization

M = as.matrix(matrix(c(2,1,1,4,1,0,-2,2,1), ncol=3, nrow=3, byrow = TRUE))
ld_factor <- function(x)
{
  counter = 2
  pos = 1
  identity = seq(1,nrow(x)*ncol(x))
  for(r in seq(nrow(x)))
  {
    for(c in seq(ncol(x)))
    {
      if(r==c)
        {identity[pos] = 1}
      else
      {identity[pos] = 0}
      pos = pos + 1
    }
  }
  
  identity = as.matrix(matrix(identity, nrow = nrow(x), ncol=ncol(x), byrow = TRUE))
  
  for (c in seq(1,ncol(x)-1)) {
    
    for( r in seq(counter, nrow(x)))
    {
      factor = -(x[r,c]/x[counter-1,c])
      new_r = factor*x[counter-1,] 
      x[r,] = new_r + x[r,]
      
      identity[r,c] = -1*factor
    }
    counter = counter + 1
  }
  
  list(x, identity)
}

ld_factor(M)

## [[1]]
##      [,1] [,2] [,3]
## [1,]    2    1    1
## [2,]    0   -1   -2
## [3,]    0    0   -4
## 
## [[2]]
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]   -1   -3    1