function  [ref_data, exp_data] = profile_cgemm_pair( ref_file, ref_model, exp_file, exp_model)
%
%   profile_sgemm_pair: compare two sgemm algorithms
%   
%  [ref_data, exp_data] = profile_cgemm_pair( ref_file, ref_model, exp_file, exp_model)
%
% Input: 
%   ref_file: input file of reference model
%   ref_name: a string described name of reference model
%   exp_file: input file of experiment model
%   exp_name: a string described name of experiment model
% 
%  format of input file, 'ref_file' and 'exp_file'
%  [dimension n]   [gpu_time (ms)]   [flops (Gflops/s)]
%
% Example: if we want to compare volkov's code and method 1, then 
%
%   profile_cgemm_pair( '../data/volkov/threads512.txt', 'volkov', '../data/method1/threads320.txt' , 'method1') 
%

str = sprintf('Load reference model %s', ref_file) ;
disp(str) ;

ref_data = load( ref_file ) ;
N = ref_data(:,1) ;
ref_gpu_time = ref_data(:,2) ;
ref_flops = ref_data(:,3) ;

str = sprintf('Load experiment model %s', exp_file) ;
disp(str) ;

exp_data = load( exp_file ) ;
exp_gpu_time = exp_data(:,2) ;
exp_flops = exp_data(:,3) ;

figure(1)
	plot(N, ref_flops, '-') ;
	hold on
		plot(N, exp_flops, '-.r') ;
	hold off
	xlabel('N') ;
	xmin = min(N) ;
	xmax = max(N) ;
	ymin = max(100, min(exp_flops)) ;
	ymax = max( [ max(exp_flops), max(ref_flops) ] ) + 50 ;
	axis([xmin xmax ymin ymax ]);
	legend( ref_model, exp_model) ;
	title('(Gflops/s)', 'fontsize', 12 ) ;
	
figure(2)
  R = exp_gpu_time./ref_gpu_time ;
	plot(N, R, '-') ;
	hold on
		plot(N, ones( length(N),1) , '-.r' ) ;
		plot(N, 0.9 *ones( length(N),1) , '-.g' ) ;
		plot(N, 0.8 *ones( length(N),1) , '-.c' ) ;		
	hold off
	xlabel('N') ;
	xmin = min(N) ;
	xmax = max(N) ;
	ymin = min(0.8, min(R)) - 0.1 ;
	ymax = max( R ) + 0.1 ;
	axis([xmin xmax ymin ymax ]);
		
	str = sprintf('R = time(%s) / time(%s)', exp_model, ref_model);
	title( str ,'fontsize', 12 ) ;	

%
% plot R(N = multiple of 8)
%
N_8 = 1024:8:max(N) ;
I_8 = extract_index( N, N_8 ) ;
figure(3)
	R_8 = R(I_8);
	plot(N(I_8), R_8, '-') ;
	hold on
		plot(N(I_8), ones( length(I_8),1) , '-.r' ) ;
	hold off
	xlabel('N') ;
	xmin = 1024 ;
	xmax = max(N) ;
	ymin = min(R_8) - 0.1 ;
	ymax = max(R_8) + 0.1 ;
	axis([xmin xmax ymin ymax ]);	
	title( 'R( N = multiple of 8)' ,'fontsize', 12 ) ;	
	
N1 = [256 512 1024 2048 4096] ;
I1 = extract_index( N, N1 ) ;	

disp(' N1 = [256 512 1024 2048 4096]') ;
ref_flops_N1 =  ref_flops( I1 ) 
exp_flops_N1 =  exp_flops( I1 ) 
R = ref_flops_N1 ./ exp_flops_N1



%
% plot R(N = multiple of 8)
%
N_64 = 1024:64:max(N) ;
I_64 = extract_index( N, N_64 ) ;
figure(4)
	plot(N(I_64), ref_flops( I_64 ), '-') ;
	hold on
		plot(N(I_64), exp_flops( I_64 ) , '-.r' ) ;
	hold off
	xlabel('N') ;
	legend( ref_model, exp_model) ;
	title( 'Gflop/s, N = multiple of 64' ,'fontsize', 12 ) ;	
