# 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. }