# Speed test of solving linear systems A*x = d

rm(list=ls()) #rm(list of objects) removes all objects from memory
graphics.off() #Closing all previously open graphs

n = 5000 #Number of linear equations and unknowns
A <- replicate(n,runif(n)) #Matrx of coefficients
d <- as.matrix(runif(n), ncol = 1) #Vector (column) of constants

# Matrix inversion
tic <- Sys.time()
x1 <- solve(A) %*% d #solve(A) is inv(A), and %*% is dot product
print(Sys.time() - tic)

# Gaussian elimination
tic <- Sys.time()
x2 <- solve(A,d)  #Left matrix division, A\d in Matlab
print(Sys.time() - tic)

X <- data.frame(x1,x2) #Comparing the two solutions
norm(x1-x2) #Differences between the two solutions