<jupyter_start><jupyter_text>Benchmarking with and without OpenBLAS, on Sparse Matrix ObjectiveCompares runtime and accuracy of computing traceinv using OpenBLAS versus not using OpenBLAS.This notebook will plot* Process time of each method.* Accuracy of computing traceinv wth respect to 128-bit data type. MatricesMatrix $\mathbf{A}$ are symmetric positive-definite matrices from real applications obtained from [sparsesuite matrix collection](https://sparse.tamu.edu/).* ['nos5'](https://sparse.tamu.edu/HB/nos5)* ['mhd4800b'](https://sparse.tamu.edu/Bai/mhd4800b)* ['bodyy6'](https://sparse.tamu.edu/Pothen/bodyy6)* ['G2_circuit'](https://sparse.tamu.edu/AMD/G2_circuit)* ['parabolic_fem'](https://sparse.tamu.edu/Wissgott/parabolic_fem)* ['StocF-1465'](https://sparse.tamu.edu/Janna/StocF-1465)* ['Bump_2911'](https://sparse.tamu.edu/Janna/Bump_2911)* ['Queen_4147'](https://sparse.tamu.edu/Janna/Queen_4147): This matrix is so large that 32-bit integers cannot hold the indices of this matrix. Rather, 32-bit `unsigned int` should be used. Consideratons**parameters:*** In SLQ and Hutchinson methods, `min_num_samples` and `max_num_samples` are fixed to 200.* In SLQ method, `lanczos_degree` is 80. Lanczos degree of 50 is no enough.* All 24 cores of Intel Xeon E5-2670 v3 processor are used for all algorithms.* The GPU results are obtained by GTX-3090 GPUs in [Genesis Cloud](https://www.genesiscloud.com/).**Notes:*** OpenBLAS only operates on dense matrices, however, here, we use sparse matrices. Here, OpenBLAS is used only on vector operations such as computing norm of a vector, or dot product. The rest of matrix-vector operations are performed without OpenBLAS. How to run the code Run on local CPURun `/benchmark/scripts/benchmark_speed.py` by    cd /benchmark/scripts    python ./benchmark_speed.py -c    where `-c` runs the code on CPU on all 32-bit, 64-bit, and 128-bit data types.The output is stored in `/benchnmark/pickle_results/benchmark_results_cpu.pickle`.The above code should be run twince,* once when `imate` is compiled using OpenBLAS (set `export USE_CBLAS=1`). Save the result into      /benchmark/pickle_results/benchmark_results_cpu_2670v3_openblas.pickle      * once when `imate` is compiled without using OpenBLAS (set `export USE_CBLAS=0`). Save the result into      /benchmark/pickle_results/benchmark_results_cpu_2670v3.pickle Run on CPU clusterTo run this script on a cluster with SLURM:    cd jobfiles    sbatch jobfile_cpu.sh   When submitting the jobs, make sure that the cpu is the same as the previous runs. For isnatnce, nodes on savio2 cluster between `n027` and `n150` are *Intel Xeon E5-2670 v3*.<jupyter_code>import numpy
import pickle
import pprint
pp = pprint.PrettyPrinter(indent=4)

from os.path import join
import matplotlib
import matplotlib.ticker
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from matplotlib.ticker import MaxNLocator
import seaborn as sns
matplotlib.rcParams['text.usetex'] = True
matplotlib.rcParams['text.latex.preamble'] = r'\usepackage{amsmath}'
#matplotlib.font_manager._rebuild() 
sns.set(font_scale=1.2)
sns.set_style("white")
sns.set_style("ticks")
plt.rc('font', family='serif')
plt.rcParams['svg.fonttype'] = 'none'
benchmark_dir = '..'
cpu_results_filename = join(benchmark_dir, 'pickle_results', 'benchmark_results_cpu_2670v3.pickle')
cpu_blas_results_filename = join(benchmark_dir, 'pickle_results',
                                 'benchmark_results_cpu_2670v3_openblas.pickle')
#cpu_blas_results_filename = join(benchmark_dir, 'pickle_results',
#                                 'benchmark_results_cpu_2680v4_openblas.pickle')


with open(cpu_results_filename, 'rb') as cpu_file:
    cpu_benchmark_result = pickle.load(cpu_file)
with open(cpu_blas_results_filename, 'rb') as cpu_blas_file:
    cpu_blas_benchmark_result = pickle.load(cpu_blas_file)
    
print('CPU Config:')
pp.pprint(cpu_benchmark_result['config'])
print('BPU BLAS Config:')
pp.pprint(cpu_blas_benchmark_result['config'])
print('')
print('CPU Devices:')
pp.pprint(cpu_benchmark_result['devices'])
print('CPU BLAS Devices:')
pp.pprint(cpu_blas_benchmark_result['devices'])

print('')
cpu_data_results = cpu_benchmark_result['data_results']
cpu_blas_data_results = cpu_blas_benchmark_result['data_results']
print('CPU num data: %d' % len(cpu_data_results))
print('CPU BLAS num data: %d' % len(cpu_blas_data_results))

cpu_nnz = numpy.zeros((len(cpu_data_results),))
cpu_size = numpy.zeros((len(cpu_data_results),))
cpu_names = []
for i in range(len(cpu_data_results)):
    cpu_nnz[i] = cpu_data_results[i]['type_results'][0]['cpu_results'][0]['info']['matrix']['nnz']
    cpu_size[i] = cpu_data_results[i]['type_results'][0]['cpu_results'][0]['info']['matrix']['size']
    cpu_name = cpu_data_results[i]['data_name']
    cpu_names.append(cpu_name.encode('unicode_escape').decode())
    
cpu_blas_nnz = numpy.zeros((len(cpu_blas_data_results),))
cpu_blas_size = numpy.zeros((len(cpu_blas_data_results),))
cpu_blas_names = []
for i in range(len(cpu_data_results)):
    cpu_blas_nnz[i] = cpu_blas_data_results[i]['type_results'][0]['cpu_results'][0]['info']['matrix']['nnz']
    cpu_blas_size[i] = cpu_blas_data_results[i]['type_results'][0]['cpu_results'][0]['info']['matrix']['size']
    cpu_blas_name = cpu_blas_data_results[i]['data_name']
    cpu_blas_names.append(cpu_blas_name.encode('unicode_escape').decode())
    

print('')
print('CPU data names:')
print(cpu_names)
print('CPU BLAS data names:')
print(cpu_blas_names)
# Sorting data by nnz, plot nnz and matrix size
sort_index = numpy.argsort(cpu_nnz)

fig, ax = plt.subplots()
x = numpy.arange(len(cpu_data_results))
ax.plot(cpu_size[sort_index], '-o', color='blue', label='size')
ax.set_ylabel('size')
ax.set_xticks(x)
#ax.set_xticklabels(names)
ax2 = ax.twinx()
ax2.plot(cpu_nnz[sort_index], '-o', color='red', label='nnz')
ax.legend(loc='upper left')
ax2.legend(loc='lower right')
ax2.set_ylabel('nnz')
ax.set_yscale('log', base=10)
ax2.set_yscale('log', base=10)

print('')
for i in range(len(cpu_data_results)):
    print('i: %d, index: %d, name: %s, nnz: %d' % (i, sort_index[i], cpu_names[sort_index[i]], cpu_nnz[i]))
def find_slope(x, y):
    """
    finds a in y = x**a
    """
    
    log2_x = numpy.log2(x)
    log2_y = numpy.log2(y)
    
    return numpy.polyfit(log2_x, log2_y, 1)[0]
## Plotting Elapsed time vs Cores

# Settings
cpu_data_result = cpu_data_results[7]
cpu_blas_data_result = cpu_blas_data_results[7]
plot_log = True
intel_color = ['#0093FF', '#0072C6', '#005EA3']         # corresponding to 32-bit, 64-bit, and 128-bit
intel_blas_color = ['#FAB300', '#C77E00']               # corresponding to 32-bit, 64-bit
markers = ['-o', '-s', '-D']                            # corresponding to 32-bit, 64-bit, and 128-bit

data_name = cpu_data_result['data_name']
print('data name: %s' % data_name)
cpu_type_results = cpu_data_result['type_results']
cpu_blas_type_results = cpu_blas_data_result['type_results']

# Normalie all wall times by cpu wall time of float type
benchmark_wall_time = cpu_data_result['type_results'][0]['cpu_results'][0]['info']['time']['alg_wall_time']

matrix_size = cpu_type_results[0]['cpu_results'][0]['info']['matrix']['size']
matrix_nnz = cpu_type_results[0]['cpu_results'][0]['info']['matrix']['nnz']
matrix_density = cpu_type_results[0]['cpu_results'][0]['info']['matrix']['density']
print('matrix size: %d' % matrix_size)
print('matrix nnz: %d' % matrix_nnz)
print('matrix density: %f' % matrix_density)

fig, ax = plt.subplots(figsize=(7, 5))

# Iterate over 32-bit, 64-bit and 128-bit data
for type_index in range(len(cpu_type_results)):
               
    cpu_type_result = cpu_type_results[type_index]
    data_type = cpu_type_result['data_type']
    cpu_results = cpu_type_result['cpu_results']

    cpu_num_threads = numpy.zeros((len(cpu_results),), dtype=int)
    cpu_tot_wall_time = numpy.zeros((len(cpu_results),), dtype=float)
    cpu_alg_wall_time = numpy.zeros((len(cpu_results),), dtype=float)
    cpu_cpu_proc_time = numpy.zeros((len(cpu_results),), dtype=float)

    for i in range(len(cpu_results)):
        cpu_num_threads[i]   = cpu_results[i]['info']['device']['num_cpu_threads']
        cpu_tot_wall_time[i] = cpu_results[i]['info']['time']['tot_wall_time']
        cpu_alg_wall_time[i] = cpu_results[i]['info']['time']['alg_wall_time']
        cpu_cpu_proc_time[i] = cpu_results[i]['info']['time']['cpu_proc_time']
        
    slope = find_slope(cpu_num_threads, cpu_alg_wall_time)
    print("CPU slope: %f" % slope)

    multi_thread = False
    if not multi_thread:
        cpu_num_cores = cpu_num_threads //2
        
    if plot_log:
        ax.plot(cpu_num_cores, cpu_alg_wall_time,
                markers[type_index], color=intel_color[type_index],
                label='no openblas ' + str(data_type) + '-bit')
    else:
        ax.plot(cpu_num_cores, benchmark_wall_time / cpu_alg_wall_time,
                markers[type_index], color=intel_color[type_index],
                label='no openblas ' + str(data_type) + '-bit')
        
    if type_index < len(cpu_blas_type_results):
        cpu_blas_type_result = cpu_blas_type_results[type_index]
        cpu_blas_results = cpu_blas_type_result['cpu_results']

        if len(cpu_blas_results) > 0:
            cpu_blas_num_threads = numpy.zeros((len(cpu_blas_results),), dtype=int)
            cpu_blas_tot_wall_time = numpy.zeros((len(cpu_blas_results),), dtype=float)
            cpu_blas_alg_wall_time = numpy.zeros((len(cpu_blas_results),), dtype=float)
            cpu_blas_cpu_proc_time = numpy.zeros((len(cpu_blas_results),), dtype=float)


            for i in range(len(cpu_blas_results)):
                cpu_blas_num_threads[i] = cpu_blas_results[i]['info']['device']['num_cpu_threads']
                cpu_blas_tot_wall_time[i]   = cpu_blas_results[i]['info']['time']['tot_wall_time']
                cpu_blas_alg_wall_time[i]   = cpu_blas_results[i]['info']['time']['alg_wall_time']
                cpu_blas_cpu_proc_time[i]   = cpu_blas_results[i]['info']['time']['cpu_proc_time']

            if not multi_thread:
                cpu_blas_num_cores = cpu_blas_num_threads //2

            if cpu_blas_num_cores.size > 1:
                slope = find_slope(cpu_blas_num_cores, cpu_blas_alg_wall_time)
                print("CPU BLAS slope: %f" % slope)

            if plot_log:
                ax.plot(cpu_blas_num_cores, cpu_blas_alg_wall_time, markers[type_index],
                    color=intel_blas_color[type_index], label='openblas' + str(data_type) + '-bit')
            else:
                ax.plot(cpu_blas_num_cores, benchmark_wall_time / cpu_blas_alg_wall_time, markers[type_index],
                    color=intel_blas_color[type_index], label='openblas ' + str(data_type) + '-bit')
        
        
# create blank rectangle
extra = Rectangle((0, 0), 1, 1, fc="w", fill=False, edgecolor='none', linewidth=0)
handles,labels = ax.get_legend_handles_labels()
handles = [extra, extra, extra, extra,
           extra, handles[0], handles[2], handles[4], 
           extra, handles[1], handles[3], extra]
labels = ["", "32 bit", "64 bit", "128 bit",
          "No BLAS", "", "", "",
          "BLAS", "", "", ""]
#handles = [handles[0], handles[2], handles[4], handles[1], handles[3]]
#labels = [labels[0], labels[2], labels[4], labels[1], labels[3]]

ax.legend(handles, labels, bbox_to_anchor=(1, 1), loc='upper left', ncol=3, handletextpad = -2)
    
if plot_log:
    plt.xscale('log', base=2)
    plt.yscale('log', base=2)
    
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
ax.get_xaxis().get_major_formatter().labelOnlyBase = False
ax.set_xticks(cpu_num_cores)

if plot_log:
    ax.set_ylabel('Elapsed time (sec)')
    ax.set_title('Elapsed Time v.s. Number of Cores')
else:
    ax.set_ylabel('Speedup')
    ax.set_title('Scalability')

ax.set_xlabel('Number of CPU cores')
ax.grid(axis='y')

# Save plot
save_dir = 'svg_plots'
filename_svg = 'benchmark_openblas_elasped_time.svg'
plt.savefig(join(benchmark_dir, save_dir, filename_svg), transparent=True, bbox_inches='tight')
print('Plot saved to "%s".' % (filename_svg))
## Plotting Elapsed time vs Matrix Size

# Settings
plot_log = True
intel_color = ['#0093FF', '#0072C6', '#005EA3']    # corresponding to 32-bit, 64-bit, and 128-bit
intel_blas_color = ['#FAB300', '#C77E00']               # corresponding to 32-bit, 64-bit
markers = ['-o', '-s', '-D']                       # corresponding to 32-bit, 64-bit, and 128-bit

# Remove some data
cpu_data_results_ = cpu_data_results.copy()
cpu_blas_data_results_ = cpu_blas_data_results.copy()

#remove_indices = [0, 1, 2]
remove_indices = []

for index in remove_indices[::-1]:
    del cpu_data_results_[index]
    del cpu_blas_data_results_[index]

# Initialize arrays
num_data = len(cpu_data_results_)
num_types = 3   # correspondng to 32-bit, 64-bit, and 128-bit
matrix_size = numpy.zeros((num_data, ), dtype=int)
matrix_name = []
matrix_nnz = numpy.zeros((num_data, ), dtype=int)
matrix_density = numpy.zeros((num_data, ), dtype=float)
data_type = numpy.zeros((num_data, num_types), dtype=int)
cpu_tot_wall_time = numpy.zeros((num_data, num_types), dtype=float)
cpu_alg_wall_time = numpy.zeros((num_data, num_types), dtype=float)
cpu_cpu_proc_time = numpy.zeros((num_data, num_types), dtype=float)
cpu_blas_tot_wall_time = numpy.zeros((num_data, num_types), dtype=float)
cpu_blas_alg_wall_time = numpy.zeros((num_data, num_types), dtype=float)
cpu_blas_cpu_proc_time = numpy.zeros((num_data, num_types), dtype=float)

for data_index in range(num_data):

    cpu_data_result = cpu_data_results_[data_index]
    matrix_name.append(cpu_data_result['data_name'])
    cpu_type_results = cpu_data_result['type_results']
    
    if data_index < len(cpu_blas_data_results_):
        cpu_blas_data_result = cpu_blas_data_results_[data_index]
        cpu_blas_type_results = cpu_blas_data_result['type_results']

    matrix_size[data_index] = cpu_type_results[0]['cpu_results'][0]['info']['matrix']['size']
    matrix_nnz[data_index] = cpu_type_results[0]['cpu_results'][0]['info']['matrix']['nnz']
    matrix_density[data_index] = cpu_type_results[0]['cpu_results'][0]['info']['matrix']['density']

    # Iterate over 32-bit, 64-bit and 128-bit data
    for type_index in range(len(cpu_type_results)):

        cpu_type_result = cpu_type_results[type_index]
        data_type[data_index, type_index] = cpu_type_result['data_type']

        cpu_results = cpu_type_result['cpu_results']

        cpu_num_threads = numpy.zeros((len(cpu_results),), dtype=int)
        cpu_tot_wall_time_ = numpy.zeros((len(cpu_results),), dtype=float)
        cpu_alg_wall_time_ = numpy.zeros((len(cpu_results),), dtype=float)
        cpu_cpu_proc_time_ = numpy.zeros((len(cpu_results),), dtype=float)

        for i in range(len(cpu_results)):
            cpu_num_threads[i]   = cpu_results[i]['info']['device']['num_cpu_threads']
            cpu_tot_wall_time_[i] = cpu_results[i]['info']['time']['tot_wall_time']
            cpu_alg_wall_time_[i] = cpu_results[i]['info']['time']['alg_wall_time']
            cpu_cpu_proc_time_[i] = cpu_results[i]['info']['time']['cpu_proc_time']
            
        cpu_max_num_threads_index = numpy.argmax(cpu_num_threads)
        #cpu_max_num_threads_index = numpy.argmin(cpu_num_threads)
        cpu_tot_wall_time[data_index, type_index] = cpu_tot_wall_time_[cpu_max_num_threads_index]
        cpu_alg_wall_time[data_index, type_index] = cpu_alg_wall_time_[cpu_max_num_threads_index]
        cpu_cpu_proc_time[data_index, type_index] = cpu_cpu_proc_time_[cpu_max_num_threads_index]
        
        if type_index < len(cpu_blas_type_results):
            cpu_blas_type_result = cpu_blas_type_results[type_index]
            cpu_blas_results = cpu_blas_type_result['cpu_results']

            if len(cpu_blas_results) > 0:
                cpu_blas_num_threads_ = numpy.zeros((len(cpu_blas_results),), dtype=int)
                cpu_blas_tot_wall_time_ = numpy.zeros((len(cpu_blas_results),), dtype=float)
                cpu_blas_alg_wall_time_ = numpy.zeros((len(cpu_blas_results),), dtype=float)
                cpu_blas_cpu_proc_time_ = numpy.zeros((len(cpu_blas_results),), dtype=float)

                for i in range(len(cpu_blas_results)):
                    cpu_blas_num_threads_[i]   = cpu_blas_results[i]['info']['device']['num_cpu_threads']
                    cpu_blas_tot_wall_time_[i] = cpu_blas_results[i]['info']['time']['tot_wall_time']
                    cpu_blas_alg_wall_time_[i] = cpu_blas_results[i]['info']['time']['alg_wall_time']
                    cpu_blas_cpu_proc_time_[i] = cpu_blas_results[i]['info']['time']['cpu_proc_time']

                cpu_max_num_cpu_threads_index = numpy.argmax(cpu_blas_num_threads_)
                cpu_blas_tot_wall_time[data_index, type_index] = \
                    cpu_blas_tot_wall_time_[cpu_max_num_cpu_threads_index]
                cpu_blas_alg_wall_time[data_index, type_index] = \
                    cpu_blas_alg_wall_time_[cpu_max_num_cpu_threads_index]
                cpu_blas_cpu_proc_time[data_index, type_index] = \
                    cpu_blas_cpu_proc_time_[cpu_max_num_cpu_threads_index]
            
    #sort_type_index = numpy.argsort(data_type[data_index, :])
    #data_type[data_index, :] = data_type[data_index, sort_type_index]
    #cpu_tot_wall_time[data_index, :] = cpu_tot_wall_time[data_index, sort_type_index]
    #cpu_alg_wall_time[data_index, :] = cpu_alg_wall_time[data_index, sort_type_index]
    #cpu_cpu_proc_time[data_index, :] = cpu_cpu_proc_time[data_index, sort_type_index]
    #cpu_blas_tot_wall_time[data_index, :] = cpu_blas_tot_wall_time[data_index, sort_type_index]
    #cpu_blas_alg_wall_time[data_index, :] = cpu_blas_alg_wall_time[data_index, sort_type_index]
    #cpu_blas_cpu_proc_time[data_index, :] = cpu_blas_cpu_proc_time[data_index, sort_type_index]
    

# Set what should be the x axis of both original and the twin plot
# All variables will be sorted by this variable in oncreasing order
#abscissa = matrix_size  
abscissa = matrix_nnz

# Set what should be on the y axis of the twin plot.
#ordinate = matrix_density
ordinate = matrix_size
#ordinate = matrix_nnz
        
sorting_index = numpy.argsort(abscissa)
abscissa = abscissa[sorting_index]
ordinate = ordinate[sorting_index]
cpu_alg_wall_time = cpu_alg_wall_time[sorting_index, :]
cpu_blas_alg_wall_time = cpu_blas_alg_wall_time[sorting_index, :]

print(matrix_name)
matrix_name_ = []
for i in sorting_index:
    matrix_name_.append(matrix_name[i])
matrix_name = matrix_name_

print(matrix_name)
print(abscissa)
            
# Plotting
fig, ax = plt.subplots(figsize=(10, 5))

#width = 0.17
width = 0.14
#logsp = numpy.logspace(-width/1.55, width/1.55, 4)
logsp = numpy.logspace(-width/1.135, width/1.135, 5)
cpu_widths = logsp[:3]
cpu_blas_widths = logsp[3:]

for type_index in range(num_types):
    
    type_ = data_type[0][type_index]
    x = cpu_widths[type_index]*abscissa
    ax.bar(x, 100.0 * cpu_alg_wall_time[:, type_index] / cpu_alg_wall_time[:, 0],
           width=width*x, edgecolor='black',
       color=intel_color[type_index], label='no BLAS ' + str(type_) + '-bit')
    
    if type_ != 128 and type_ != 0:
        x = cpu_blas_widths[type_index]*abscissa
        ax.bar(x, 100.0 * cpu_blas_alg_wall_time[:, type_index] / cpu_alg_wall_time[:, 0],
               width=width*x, edgecolor='black',
           color=intel_blas_color[type_index], label='BLAS ' + str(type_) + '-bit')

"""
for type_index in range(num_types):
   
    type_ = data_type[0, type_index]
    ax.plot(abscissa, cpu_alg_wall_time[:, type_index],
                markers[type_index], color=intel_color[type_index],
            label='no BLAS ' + str(type_) + '-bit')

    if type_ != 128:
        ax.plot(abscissa, cpu_blas_alg_wall_time[:, type_index], markers[type_index],
                color=intel_blas_color[type_index],
                label='BLAS ' + str(type_) + '-bit')
"""
        
# create blank rectangle
extra = Rectangle((0, 0), 1, 1, fc="w", fill=False, edgecolor='none', linewidth=0)
handles,labels = ax.get_legend_handles_labels()
handles = [extra, extra, extra, extra,
           extra, handles[0], handles[2], handles[4], 
           extra, handles[1], handles[3], extra]
labels = ["", "32 bit", "64 bit", "128 bit",
          "no BLAS", "", "", "",
          "BLAS", "", "", ""]


ax2 = ax.twinx()
ax2.plot(abscissa, ordinate, 'o', color='gray', markerfacecolor='white', label='matrix size')
#ax.legend(handles, labels, bbox_to_anchor=(1.1, 1), loc='upper left', ncol=3, handletextpad = -2,
#          fontsize='small')
ax.legend(handles, labels, loc='upper left', ncol=3, handletextpad = -2, fontsize='small')
#ax.legend(handles, labels, bbox_to_anchor=(1.02, 1), loc='upper left', ncol=3, handletextpad = -2)
#ax2.legend(bbox_to_anchor=(1.15, 0.7), loc='upper left', ncol=1, fontsize='small')
ax2.legend(bbox_to_anchor=(0, 0.7), loc='upper left', ncol=1, fontsize='small')
    
plt.xscale('log', base=10)
#ax.set_yscale('log', base=10)
ax2.set_yscale('log', base=10)
ax2.set_ylim([1e2, 1e7])
ax.yaxis.set_major_formatter(matplotlib.ticker.PercentFormatter(decimals=0, xmax=100))

#ax.xaxis.set_major_locator(MaxNLocator(integer=True))
#ax.get_xaxis().get_major_formatter().labelOnlyBase = False
#ax.set_xticks(cpu_blas_num_cpu_devices)

ax.set_ylabel('Elapsed time (relative to 32-bit no BLAS)')
ax2.set_ylabel('Matrix size')
ax.set_title('Elapsed Time v.s. Sparse Matrix Size')
#ax.set_xlabel('Matrix size')
ax.set_xlabel('Matrix nnz')
ax.grid(axis='y')

# Save plot
save_dir = 'svg_plots'
filename_svg = 'benchmark_openblas_sparse_elasped_time.svg'
plt.savefig(join(benchmark_dir, save_dir, filename_svg), transparent=True, bbox_inches='tight')
print('Plot saved to "%s".' % (filename_svg))
## Plotting Elapsed time vs Matrix Size
# Note: The errors are computed by comparing with 128-bit with no CBLAS.

# Settings
plot_log = True
intel_color = ['#0093FF', '#0072C6', '#005EA3']    # corresponding to 32-bit, 64-bit, and 128-bit
intel_blas_color = ['#FAB300', '#C77E00']          # corresponding to 32-bit, 64-bit
markers = ['-o', '-s', '-D']                       # corresponding to 32-bit, 64-bit, and 128-bit

# Remove some data
cpu_data_results_ = cpu_data_results.copy()
cpu_blas_data_results_ = cpu_blas_data_results.copy()

remove_indices = [0, 1, 2]
#remove_indices = []

for index in remove_indices[::-1]:
    del cpu_data_results_[index]
    del cpu_blas_data_results_[index]

# Get number of cpu and cpu blas results from one of the datasets
num_cpu_results = len(cpu_data_results_[0]['type_results'][0]['cpu_results'])
num_cpu_blas_results = len(cpu_blas_data_results_[0]['type_results'][0]['cpu_results'])


# Initialize arrays
num_data = len(cpu_data_results_)
num_types = 3   # correspondng to 32-bit, 64-bit, and 128-bit
matrix_size = numpy.zeros((num_data, ), dtype=int)
matrix_name = []
matrix_nnz = numpy.zeros((num_data, ), dtype=int)
matrix_density = numpy.zeros((num_data, ), dtype=float)
data_type = numpy.zeros((num_types, ), dtype=int)
cpu_trace = numpy.zeros((num_data, num_types, num_cpu_results), dtype=float)
cpu_error = numpy.zeros((num_data, num_types, num_cpu_results), dtype=float)
cpu_blas_trace = numpy.zeros((num_data, num_types, num_cpu_blas_results), dtype=float)
cpu_blas_error = numpy.zeros((num_data, num_types, num_cpu_blas_results), dtype=float)

cpu_blas_trace[:] = numpy.nan
cpu_blas_error[:] = numpy.nan


for data_index in range(num_data):

    cpu_data_result = cpu_data_results_[data_index]
    matrix_name.append(cpu_data_result['data_name'])
    cpu_type_results = cpu_data_result['type_results']
    
    if data_index < len(cpu_blas_data_results):
        cpu_blas_data_result = cpu_blas_data_results_[data_index]
        cpu_blas_type_results = cpu_blas_data_result['type_results']

    matrix_size[data_index] = cpu_type_results[0]['cpu_results'][0]['info']['matrix']['size']
    matrix_nnz[data_index] = cpu_type_results[0]['cpu_results'][0]['info']['matrix']['nnz']
    matrix_density[data_index] = cpu_type_results[0]['cpu_results'][0]['info']['matrix']['density']

    # Iterate over 32-bit, 64-bit and 128-bit data
    for type_index in range(len(cpu_type_results)):

        cpu_type_result = cpu_type_results[type_index]
        data_type[type_index] = cpu_type_result['data_type']
        cpu_results = cpu_type_result['cpu_results']

        cpu_num_threads_ = numpy.zeros((len(cpu_results),), dtype=int)

        for i in range(len(cpu_results)):
            cpu_num_threads[i] = cpu_results[i]['info']['device']['num_cpu_threads']
            cpu_trace[data_index, type_index, i] = cpu_results[i]['trace']
            cpu_error[data_index, type_index, i] = cpu_results[i]['info']['error']['absolute_error']
            
        cpu_sorting_index = numpy.argsort(cpu_num_threads)
        cpu_trace[:, :, :] = cpu_trace[:, :, cpu_sorting_index]
        cpu_error[:, :, :] = cpu_error[:, :, cpu_sorting_index]
        
        if type_index < len(cpu_blas_type_results):
            cpu_blas_type_result = cpu_blas_type_results[type_index]
            cpu_blas_results = cpu_blas_type_result['cpu_results']

            if len(cpu_blas_results) > 0:
                cpu_blas_num_threads = numpy.zeros((len(cpu_blas_results),), dtype=int)

                for i in range(len(cpu_blas_results)):
                    cpu_blas_num_threads[i] = cpu_blas_results[i]['info']['device']['num_cpu_threads']
                    cpu_blas_trace[data_index, type_index, i] = cpu_blas_results[i]['trace']
                    cpu_blas_error[data_index, type_index, i] = \
                        cpu_blas_results[i]['info']['error']['absolute_error']

                cpu_blas_sorting_index = numpy.argsort(cpu_blas_num_threads)
                cpu_blas_trace[:, :, :] = cpu_blas_trace[:, :, cpu_blas_sorting_index]
                cpu_blas_error[:, :, :] = cpu_blas_error[:, :, cpu_blas_sorting_index]

cpu_trace = numpy.mean(cpu_trace, axis=2)
cpu_error = numpy.mean(cpu_error, axis=2)
cpu_blas_trace = numpy.mean(cpu_blas_trace, axis=2)
cpu_blas_error = numpy.mean(cpu_blas_error, axis=2)

# Relative comparison w.r.t 128-bit (in percent)
cpu_trace_relative = numpy.zeros((num_data, num_types), dtype=float)
cpu_blas_trace_relative = numpy.zeros((num_data, num_types), dtype=float)
cpu_error_relative = numpy.zeros((num_data, num_types), dtype=float)
cpu_blas_error_relative = numpy.zeros((num_data, num_types), dtype=float)
percent = 100.0


for data_index in range(num_data):
    benchmark_cpu_trace = numpy.abs(cpu_trace[data_index, -1])
    benchmark_cpu_blas_trace = numpy.abs(cpu_blas_trace[data_index, -1])
    
    for type_index in range(num_types):
        cpu_trace_relative[data_index, type_index ] = percent * \
            numpy.abs(cpu_trace[data_index, type_index] - benchmark_cpu_trace) / benchmark_cpu_trace
        cpu_blas_trace_relative[data_index, type_index ] = percent * \
            numpy.abs(cpu_blas_trace[data_index, type_index] - benchmark_cpu_trace) / benchmark_cpu_trace
        
        cpu_error_relative[data_index, type_index ] = percent * \
            numpy.abs(cpu_error[data_index, type_index]) / benchmark_cpu_trace
        cpu_blas_error_relative[data_index, type_index ] = percent * \
            numpy.abs(cpu_blas_error[data_index, type_index]) / benchmark_cpu_trace

    
# Plotting
fig, ax = plt.subplots(figsize=(7, 5))

# Set what should be the x axis of both original and the twin plot
# All variables will be sorted by this variable in oncreasing order
#abscissa = matrix_size  
abscissa = matrix_nnz

width = 0.17
logsp = numpy.logspace(-width/1.55, width/1.55, 4)
cpu_widths = logsp[:2]
cpu_blas_widths = logsp[2:]

for type_index in range(num_types-1):
    
    type_ = data_type[type_index]
    x = cpu_widths[type_index]*abscissa
    ax.bar(x, cpu_trace_relative[:, type_index], yerr=cpu_error_relative[:, type_index],
           width=width*x, edgecolor='black',
       color=intel_color[type_index], label='no BLAS ' + str(type_) + '-bit')
    
    if type_ != 128:
        x = cpu_blas_widths[type_index]*abscissa
        ax.bar(x, cpu_blas_trace_relative[:, type_index], yerr=cpu_blas_error_relative[:, type_index],
               width=width*x, edgecolor='black',
           color=intel_blas_color[type_index], label='BLAS ' + str(type_) + '-bit')

        
# create blank rectangle
extra = Rectangle((0, 0), 1, 1, fc="w", fill=False, edgecolor='none', linewidth=0)
handles,labels = ax.get_legend_handles_labels()
handles = [extra, extra, extra,
           extra, handles[0], handles[2], 
           extra, handles[1], handles[3]]
labels = ["", "32 bit", "64 bit",
          "no BLAS", "", "",
          "BLAS", "", ""]
#handles = [handles[0], handles[2], handles[4], handles[1], handles[3]]
#labels = [labels[0], labels[2], labels[4], labels[1], labels[3]]

#ax.legend(handles, labels, bbox_to_anchor=(1, 1), loc='upper left', ncol=3, handletextpad = -2,
#          fontsize='small')
ax.legend(handles, labels, loc='upper left', ncol=3, handletextpad = -2, fontsize='small')
    
plt.xscale('log', base=10)
plt.yscale('log', base=10)

#ax.xaxis.set_major_locator(MaxNLocator(integer=True))
#ax.get_xaxis().get_major_formatter().labelOnlyBase = False
#ax.set_xticks(data_type)
#ax.yaxis.set_major_formatter(matplotlib.ticker.PercentFormatter(decimals=0, xmax=100))
ax.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(lambda x, pos: '$10^{%d}\%%$' % numpy.log10(x)))
ax.set_ylim([1e-2, 1e2])
#ax.set_ylim([0, 30])


ax.set_ylabel('Relative error compared to 128-bit precision')
ax.set_title('Accuracy of estimating $\mathrm{trace}(\mathbf{A}^{-1})$')

ax.set_xlabel('Matrix nnz')
ax.grid(axis='y')

# Save plot
save_dir = 'svg_plots'
filename_svg = 'benchmark_openblas_spaese_accuracy.svg'
plt.savefig(join(benchmark_dir, save_dir, filename_svg), transparent=True, bbox_inches='tight')
print('Plot saved to "%s".' % (filename_svg))
data_index = 3

print('32-bit')
print(cpu_data_results[data_index]['type_results'][0]['cpu_results'][0]['trace'])
print(cpu_data_results[data_index]['type_results'][0]['cpu_results'][1]['trace'])
print(cpu_data_results[data_index]['type_results'][0]['cpu_results'][2]['trace'])
print('\n64-bit')
print(cpu_data_results[data_index]['type_results'][1]['cpu_results'][0]['trace'])
print(cpu_data_results[data_index]['type_results'][1]['cpu_results'][1]['trace'])
print(cpu_data_results[data_index]['type_results'][1]['cpu_results'][2]['trace'])
print('\n128-bit')
print(cpu_data_results[data_index]['type_results'][2]['cpu_results'][0]['trace'])
print(cpu_data_results[data_index]['type_results'][2]['cpu_results'][1]['trace'])
print(cpu_data_results[data_index]['type_results'][2]['cpu_results'][2]['trace'])
data_index = 3

print('32-bit')
print(cpu_blas_data_results[data_index]['type_results'][0]['cpu_results'][0]['trace'])
print(cpu_blas_data_results[data_index]['type_results'][0]['cpu_results'][1]['trace'])
print(cpu_blas_data_results[data_index]['type_results'][0]['cpu_results'][2]['trace'])
print('\n64-bit')
print(cpu_blas_data_results[data_index]['type_results'][1]['cpu_results'][0]['trace'])
print(cpu_blas_data_results[data_index]['type_results'][1]['cpu_results'][1]['trace'])
print(cpu_blas_data_results[data_index]['type_results'][1]['cpu_results'][2]['trace'])<jupyter_output>32-bit
20828.322
20865.422
20821.072

64-bit
20582.47750182964
20594.867831432584
20578.04029093108
{'fraction_non_alphanumeric': 0.1419311791059003, 'fraction_numerical': 0.021605820862326437, 'mean_word_length': 4.769752520623282, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 10, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '391555', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 11728, 'n_tokens_neox': 11283, 'n_words': 1913, 'path': 'benchmark/notebooks/plot_benchmark_openblas_sparse.ipynb', 'repo_name': 'ameli/TraceInv', 'size': '31473'}