# Measuring the floating point performance of your machine

Matrix-Matrix multiplication is often used to benchmark machines because the mathematics is such that it is one of the few operations where one can obtain close to theoretical peak performance in pratice.

The number of floating point operations(Flops) in a Matrix-Matrix multiplication of two $N$ x $N$ matrices is given by $2N^3-N^2$ (See http://www2.hawaii.edu/~norbert/CompPhys/compphysicsch10.html for details of how this was determined).

For this benchmark, we construct two random $N$ x $N$ matrices and time how long it takes to multiply them together, $t$. We then compute the number of Gigaflops/s exhibited by that operation via $\frac{2N^3-N^2}{t*10^9}$. We do this several times for any given $N$ and record the average and the best (peak) results obtained. $N$ is then varied from 250 to 10,000 in steps of 250 to see how performance varies with matrix size. Usually, the highest performance is obtained by multiplying large matrices together.

For highest performance, you should use a version of numpy that has been linked against a high performance BLAS library such as OpenBLAS or the Intel MKL(https://software.intel.com/en-us/intel-mkl). The Anaconda Python distribution includes the Intel MKL by default on Windows and Linux (Mac includes its own high performance BLAS library).

In [1]:
import timeit
import matplotlib.pyplot as plt

def bench_mmm(startn,maxn,step,loops):
count=0

#Preallocate results lists
avg_gflops = int((1+(maxn-startn)/step))*[0]
peak_gflops = int((1+(maxn-startn)/step))*[0]
raw_times = [int(loops)*[0] for i in range(int(1+(maxn-startn)/step))]
all_gflops = [int(loops)*[0] for i in  range(int(1+(maxn-startn)/step))]
mat_size=int((1+(maxn-startn)/step))*[0]

for n in range(startn,maxn+step,step):
setup_string = "from pylab import rand,dot;n=%d;a=rand(n,n);b=rand(n,n)" % n
time_list = timeit.repeat("a.dot(b)", setup=setup_string, repeat=loops,number=1)
raw_times[count] = time_list
total_time = sum(time_list)
avg_time = total_time / loops
peak_time = min(time_list)
num_ops = 2*n**3-n**2
avg_gflops[count] = (num_ops/avg_time)/10**9
peak_gflops[count] = (num_ops/peak_time)/10**9
all_gflops[count] = [(num_ops/time)/10**9 for time in raw_times[count]]
mat_size[count] = n
count=count+1

plt.plot(mat_size,avg_gflops,'*-',label="Average over %d runs" %loops)
plt.plot(mat_size,peak_gflops,'*-',label="Peak")
plt.xlabel('Matrix Size');
plt.ylabel('GFlop/s');
plt.show()

return(max(peak_gflops),raw_times,all_gflops)

peak_flops = bench_mmm(250,10000,250,5)

--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) <ipython-input-1-428b823d31f7> in <module>() 35 return(max(peak_gflops),raw_times,all_gflops) 36 ---> 37 peak_flops = bench_mmm(250,10000,250,5) <ipython-input-1-428b823d31f7> in bench_mmm(startn, maxn, step, loops) 14 for n in range(startn,maxn+step,step): 15 setup_string = "from pylab import rand,dot;n=%d;a=rand(n,n);b=rand(n,n)" % n ---> 16 time_list = timeit.repeat("a.dot(b)", setup=setup_string, repeat=loops,number=1) 17 raw_times[count] = time_list 18 total_time = sum(time_list) /usr/lib/python3.4/timeit.py in repeat(stmt, setup, timer, repeat, number) 224 repeat=default_repeat, number=default_number): 225 """Convenience function to create Timer object and call repeat method.""" --> 226 return Timer(stmt, setup, timer).repeat(repeat, number) 227 228 def main(args=None, *, _wrap_timer=None): /usr/lib/python3.4/timeit.py in repeat(self, repeat, number) 212 r = [] 213 for i in range(repeat): --> 214 t = self.timeit(number) 215 r.append(t) 216 return r /usr/lib/python3.4/timeit.py in timeit(self, number) 184 gc.disable() 185 try: --> 186 timing = self.inner(it, self.timer) 187 finally: 188 if gcold: /usr/lib/python3.4/timeit.py in inner(_it, _timer) KeyboardInterrupt: 
In [ ]:
#Maximum flops found
peak_flops[0]

In [ ]:
#Run it again
peak_flops = bench_mmm(250,10000,250,5)

In [ ]:
peakflops