using LinearAlgebra

using BenchmarkTools # a useful Julia package for performance benchmarking

n = [10,100,500,1000,2000]
LinearAlgebra.BLAS.set_num_threads(1) # benchmarking on multiple cores is weird
t = [@belapsed(lu($(rand(n,n))), evals=1) for n in n]

using PyPlot
loglog(n, t*1e9, "bo-")
loglog(n, n.^3, "k--")
xlabel("matrix size n")
ylabel("time (ns)")
legend(["time", L"n^3"])
title("time for LU factorization")

ts = [@belapsed($(lu(rand(n,n))) \ $(rand(n))) for n in n]

loglog(n, ts*1e9, "bo-")
loglog(n, n.^2, "k--")
xlabel("matrix size n")
ylabel("time (ns)")
legend(["time", L"n^2"])
title("time for LU solve")

t[end], n[end] # the last measured time and n for LU factorization

secs = t[end] * (1e6/n[end])^3 # this many seconds

# convert to a human time period
using Dates
Dates.canonicalize(Dates.CompoundPeriod(Dates.Second(round(Int,secs))))

println((1e6)^2 * sizeof(Float64) / 2^30, " GiB for a 10⁶×10⁶ matrix")