Advertisement
DrNostril

Untitled

Sep 5th, 2021 (edited)
216
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # input a square matrix!
  2. LUdecomposition <- function(m) {
  3.  
  4.   # throw an error if input is not square
  5.   dimensions <- dim(m)
  6.   if (dimensions[1] != dimensions[2]) {
  7.     return("Matrix not square!")
  8.  
  9.   } else {
  10.   # calculate lower and upper triangles
  11.     d <- dimensions[1]
  12.     L <- diag(d)
  13.     U <- m
  14.  
  15.     # iterate over rows 2 through d
  16.     for (x in 1:d) {
  17.  
  18.       # make sure m[1,1] = 1
  19.       if (U[1,1] != 1) {
  20.         U[1,] <- U[1,] / U[1,1]
  21.       }
  22.      
  23.       if (x > 1) {
  24.       # continue decomposing the upper triangle
  25.         for (y in 1:(x-1)) {
  26.          
  27.           # avoid dividing by zero
  28.           nonzero <- 1
  29.           for (z in (x-1):1) {
  30.             if (U[z,y] != 0) {
  31.               nonzero <- z
  32.               break
  33.             }
  34.           }
  35.          
  36.           alpha = U[x,y] / U[nonzero,y]
  37.           U[x,] <- U[x,] - (alpha * U[nonzero,])
  38.        
  39.         # add the multiplier to the lower triangle
  40.           L[x,y] <- alpha
  41.         }
  42.       }
  43.     }
  44.   }
  45.     return(list(L, U))
  46. }
Advertisement
RAW Paste Data Copied
Advertisement