###############################################################################
#
# LinpackR
# This file contains a version of the Linpack benchmark program written in the
# R programming language.
# More details are available here:  z.umn.edu/linpackr
#
# Copyright 2017 David J. Lilja
# Contact:  www.arctic.umn.edu
# This work is licensed under the Creative Commons Attribution-NonCommercial 4.0
# International License.  
# You are free to:
# Share - copy and redistribute the material in any medium or format
# Adapt - remix, transform, and build upon the material
# Under the following terms:
# Attribution – You must give appropriate credit, provide a link to the license, and 
# indicate if changes were made. You may do so in any reasonable manner, but not in any 
# way that suggests the licensor endorses you or your use.
# NonCommercial – You may not use the material for commercial purposes.
# 
# Attribution:
# See z.umn.edu/linpackr for details on how to appropriately cite this work
# for attribution.
# 
# Although every precaution has been taken to verify the accuracy of the information 
# contained herein, the author and publisher assume no responsibility for any errors or 
# omissions. No  liability is assumed for damages that may result from the use of information 
# contained within.
# 
###############################################################################


###############################################################################
######## Start of the function run_solver() ###################################
###############################################################################
# This function does the work of solving the linear systems of equations
# and calculating the performance metric value.

run_solver <- function(n) {

	# Generate the A matrix and b vector with random values between -1 and +1.
	# n is the size of the system of equations to be solved.
	A <- matrix(runif(n*n, min=-1, max=1) ,n)
	b <- matrix(runif(n, min=-1, max = 1) ,n)

	# Measure the time required to solve the system Ax = b for x.
	exec_time <- system.time(
	x <- solve(A,b),
        gcFirst=TRUE
	)[3]

	# Compute the number of flops using the standard Linpack formula
	flops <- (2/3)*n^3 + 2*n^2

	# Compute GFLOPS
	gflops <- (flops / exec_time) / 1e9

	# Compute the numerical validation of the solution - it should be < O(1)
	validation <- norm(A %*% x - b) / (norm(A) * norm(x) * n * .Machine$double.eps)

	# Print the results for this value of n
	cat(sprintf("N = %5d,",n))
	cat(sprintf(" GFLOPS = %5.1f,",gflops))
	cat(sprintf(" Validation = %6.4f,",validation))
	cat(sprintf(" Total execution time = %5.1f",exec_time))
	cat(sprintf("\n"))
}
###############################################################################
######## End of the function run_solver() #####################################
###############################################################################


###############################################################################
######## Start of the main body of the program ################################
###############################################################################
# This is the main body of the program.  You can set a few parameters here to change
# the size of the problem, the number of threads, and to print the results to a file.

# Change the following line to the list of problem sizes you want to measure.
problem_sizes <- c(1000, 2000, 3000, 4000, 5000, 10000)

# If you are running the Microsoft R Open version of the R environment
# (https://mran.revolutionanalytics.com/rro/)
# you can change the number of threads used to execute the program.  Uncomment
# the following line to set the number of threads you want.  If you do nothing,
# the system will typically use the maximum number of threads available.
# This example sets the number of threads to two:
# setMKLthreads(2)

# Uncomment the following line to direct the output to the desired file
# and simultaneously to the console.
sink("linpackR-output.txt", split=TRUE)

# Now run the solver for each value of n selected above.
for (n in problem_sizes) {
	run_solver(n)
}

# Uncomment this line if you uncommented the sink(...) function above.
sink()

###############################################################################
######## End of the main body of the program ##################################
###############################################################################