% specs for K40 GPU (note: K80s are 2x K40s)
bw = 240; % GB/s in single precision
flops = 4.373 * 1e3/2; % gigaflops in single precision

% axes of roofline plot
roofY = [0, flops, flops]; 
roofX = [0, flops/bw, 60]; % max (GFLOPS/s) / (GB/s) bandwidth 

%% measured GFLOPS/bw for original and final optimized reduction kernel

microsec_runtime = [273.91 48.862]; % microseconds
bw = [19.127+193.29/1000 105.22+431.94/1000]; 
sec_runtime = microsec_runtime/1e6; 
gflops_reduc = ([1040384 1572864]/1e9)./sec_runtime; 
arith_intensity_reduc = gflops_reduc./bw;

%% axpy kernels 

% axpy kernels with 1 iteration (chunk size 8)
bw = [28.239+57.336, 50.676+102.76];
microsec_runtime = [200 95];
sec_runtime = microsec_runtime/1e6; 
gflops_axpy1 = ([2097152 2097152]/1e9)./sec_runtime; 
arith_intensity_axpy1 = gflops_axpy1./bw;

% axpy kernels with 10 iterations (chunk size 8)
bw = [75, 150];
microsec_runtime = [226.14, 98.493];
sec_runtime = microsec_runtime/1e6; 
gflops_axpy2 = ([11534336 11534336]/1e9)./sec_runtime; 
arith_intensity_axpy2 = gflops_axpy2./bw;

% axpy kernels with 50 iterations (chunk size 8)
bw = [33.827+16.857, 81.415+40.810];
microsec_runtime = [324.98, 123.61];
sec_runtime = microsec_runtime/1e6; 
gflops_axpy3 = ([53477376 53477376]/1e9)./sec_runtime; 
arith_intensity_axpy3 = gflops_axpy3./bw;

%% measured GFLOPS/bw for naive matrix-matrix multiplication kernel (N=1024)

bw = [211.66/1000+7.5254, 11.691+20.311, 1.7706+28.749];
millisec_runtime = [24.694 1.2593 439.60/1e3];
sec_runtime = millisec_runtime/1e3;
gflops_matmat = ([2148532224 2016215040 134414336]/1e9)./sec_runtime;
arith_intensity_matmat = gflops_matmat./bw;

%% plot points on roofline plot

plot(roofX,roofY,'linewidth',2)
hold on
plot(arith_intensity_reduc,gflops_reduc,'.','linewidth',2,'markersize',12)
text(arith_intensity_reduc,gflops_reduc,{'Naive reduce','Opt reduce'})

plot(arith_intensity_axpy1,gflops_axpy1,'^','linewidth',2,'markersize',12)
text(arith_intensity_axpy1,gflops_axpy1,{'axpy (1 iter)','Opt axpy (1 iter)'})
plot(arith_intensity_axpy2,gflops_axpy2,'s','linewidth',2,'markersize',12)
text(arith_intensity_axpy2,gflops_axpy2,{'axpy (10 iter)','Opt axpy (10 iter)'})
plot(arith_intensity_axpy3,gflops_axpy3,'x','linewidth',2,'markersize',12)
text(arith_intensity_axpy3,gflops_axpy3,{'axpy (50 iter)','Opt axpy (50 iter)'})

plot(arith_intensity_matmat,gflops_matmat,'o','linewidth',2,'markersize',12)
text(arith_intensity_matmat,gflops_matmat,{'Naive matmult','sgemm 1','sgemm 2'})

set(gca,'fontsize',14)
xlabel('Arithmetic intensity','fontsize',14)
ylabel('GFLOPS/sec','fontsize',14)
grid on
