HPC/Project1/project_1_maggioni_claudio.tex

\documentclass[unicode,11pt,a4paper,oneside,numbers=endperiod,openany]{scrartcl}

\input{assignment.sty}
\usepackage{float}
\usepackage{subcaption}
\usepackage{graphicx}
\usepackage{tikz}
\usetikzlibrary{decorations.markings}
\usepackage{url}
\hypersetup{pdfborder = {0 0 0}}
\usepackage{xcolor}
\usepackage{multirow}
\usepackage{makecell}
\usepackage{booktabs}
\usepackage{algorithm2e}
\usepackage[nomessages]{fp}

\begin{document}

\setassignment
\setduedate{12.10.2022 (midnight)}

\serieheader{High-Performance Computing Lab}{2022}{Student: Claudio Maggioni}{
    Discussed with: M. Cattaneo, G. De Vita, V. Karpenko}{Solution for Project 1}{}
\newline

%\assignmentpolicy
%In this project you will practice memory access optimization,
%performance-oriented programming, and OpenMP parallelizaton on the ICS Cluster.  

\tableofcontents

\section{Explaining Memory Hierarchies \punkte{25}}

\subsection{Memory Hierarchy Parameters of the Cluster}
By invoking \texttt{likwid-topology} for the cache topology and \texttt{free -g}
for the amount of primary memory, the following memory hierarchy parameters are
found:

\begin{center}
\begin{tabular}{llll}
Main memory & 62 GB              \\
L3 cache    & 25 MB per socket   \\
L2 cache    & 256 kB per core    \\
L1 cache    & 32 kB per core    
\end{tabular}
\end{center}

All values are reported using base 2 IEC byte units. The cluster has 2 sockets
and a total of 20 cores (10 per socket). The cache topology diagram reported by
\texttt{likwid-topology -g} is shown in Figure \ref{fig:topo}.

\pagebreak[4]

\begin{figure}[t]
    \begin{center}
       Socket 0:\vspace{0.3cm}
        \begin{tabular}{|l|l|l|l|l|l|l|l|l|l|}
            \hline 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\\hline
            32 kB & 32 kB & 32 kB & 32 kB & 32 kB & 32 kB & 32 kB & 32 kB & 32
            kB & 32 kB \\\hline
            256 kB & 256 kB & 256 kB & 256 kB & 256 kB & 256 kB & 256 kB & 256
            kB & 256 kB & 256 kB \\\hline
            \multicolumn{10}{|c|}{25 MB} \\\hline
      \end{tabular}\vspace{0.8cm}\\
       Socket 1:\vspace{0.3cm}
       \begin{tabular}{|l|l|l|l|l|l|l|l|l|l|}
            \hline 10 & 11 & 12 & 13 & 14 & 15 & 16 & 17 & 18 & 19 \\\hline
            32 kB & 32 kB & 32 kB & 32 kB & 32 kB & 32 kB & 32 kB & 32 kB & 32
            kB & 32 kB \\\hline
            256 kB & 256 kB & 256 kB & 256 kB & 256 kB & 256 kB & 256 kB & 256
            kB & 256 kB & 256 kB \\\hline
            \multicolumn{10}{|c|}{25 MB} \\\hline
       \end{tabular}
    \end{center}
    \caption{Cache topology diagram as outputted by \texttt{likwid-topology -g}.
    Byte sizes all in IEC units.}
    \label{fig:topo}
\end{figure}

\subsection{Memory Access Pattern of \texttt{membench.c}}

\begin{figure}[t]
\begin{center}
\begin{tikzpicture}
    \tikzset{->-/.style={decoration={
  markings,
  mark=at position .75 with {\arrow{>}}},postaction={decorate}}};

    \draw (0,0) grid (5,1);
    \draw [dashed] (5,0) -- (5.5,0);
    \draw [dashed] (5,1) -- (5.5,1);
    \draw [dashed] (6.5,0) -- (7,0);
    \draw [dashed] (6.5,1) -- (7,1);
    \draw (7,0) grid (12,1);

    \foreach \r in {0,1,...,4}{
        \fill (\r + 0.5,0.5) circle [radius=2pt];
        \draw[->-] (\r-0.5,0.5) to[bend left] (\r+0.5,0.5);
        \draw (\r + 0.5, -0.5) node {$\r$};
    }
    \draw[->-] (4.5,0.5) to[bend left] (5.5,0.5);
    \foreach \r in {7,8,...,11}{
        \fill (\r + 0.5,0.5) circle [radius=2pt];
        \FPeval{l}{round(\r + 128 - 12, 0)}
        \draw[->-] (\r-0.5,0.5) to[bend left] (\r+0.5,0.5);
        \draw (\r + 0.5, -0.5) node {$\l$};
    }
    
    \draw (0,-3) grid (3,-2);
    \draw [dashed] (3,-2) -- (3.5,-2);
    \draw [dashed] (3,-3) -- (3.5,-3);
    \draw [dashed] (4,-2) -- (4.5,-2);
    \draw [dashed] (4,-3) -- (4.5,-3);
    \draw (4.5,-2) -- (7.5,-2);
    \draw (4.5,-3) -- (7.5,-3);
    \foreach \r in {4.5,5.5,...,7.5}{
        \draw (\r,-3) -- (\r,-2);
    }
    \draw [dashed] (7.5,-2) -- (8,-2);
    \draw [dashed] (7.5,-3) -- (8,-3);
    \draw [dashed] (8.5,-2) -- (9,-2);
    \draw [dashed] (8.5,-3) -- (9,-3);
    \draw (9,-3) grid (12,-2);

    \fill (0.5,-2.5) circle [radius=2pt];
    \fill (6,-2.5) circle [radius=2pt];
    \foreach \r in {0,1,2}{
        \draw (\r + 0.5, -3.5) node {$\r$};
    }
    \foreach \r in {9,10,11}{
        \FPeval{l}{round(\r - 12, 0)}
        \draw (\r + 0.5, -3.5) node {\tiny $2^{20}  \l$};
    }
    \draw (5, -3.5) node {\tiny $2^{19} - 1$};
    \draw (6, -3.5) node {\tiny $2^{19}$};
    \draw (7,-3.5) node {\tiny $2^{19} + 1$};
    \draw[->-] (-0.5,-2.5) to[bend left] (0.5,-2.5);
    \draw[->-] (0.5,-2.5) to[bend left] (6,-2.5);
\end{tikzpicture}
\end{center}
    \caption{Memory access patterns of \texttt{membench.c} for \texttt{csize =
    128} and \texttt{stride = 1} (above) and for \texttt{csize = $2^{20}$} and
    \texttt{stride = $2^{19}$} (below)}
    \label{fig:access}
\end{figure}

The benchmark \texttt{membench.c} measures the average time of repeated read and
write overations across a set of indices of a stack-allocated array of 32-bit
signed integers. The indices vary according to the access pattern used, which in
turn is defined by two variables, \texttt{csize} and \texttt{stride}.
\texttt{csize} is an upper bound on the index value, i.e.  (one more of) the
highest index used to access the array in the pattern.  \texttt{stride}
determines the difference between array indexes over access iterations, i.e. a
\texttt{stride} of 1 will access every array index, a \texttt{stride} of 2 will
skip every other index, a \texttt{stride} of 4 will access one index then skip 3
and so on and so forth. The benchmark stops when the index to access is strictly
greater than \texttt{csize - stride}.

Therefore, for \texttt{csize = 128} and \texttt{stride = 1} the array will
access all indexes between 0 and 127 sequentially, and for \texttt{csize =
$2^{20}$} and \texttt{stride = $2^{19}$} the benchmark will access index 0, then
index $2^{19}-1$. The access patterns for these
two configurations are shown visually in Figure \ref{fig:access}.

By running the \texttt{membench.c} both on my personal laptop and on the
cluster, the results shown in Figure \ref{fig:mem} are obtained. \textit{csize}
values are shown as different data series and labeled by byte size and
\textit{stride} values are mapped on the $x$ axis by the byte-equivalent value
as well\footnote{Byte values are a factor of 4 greater than the values used in
the code and in Figure \ref{fig:mem}. This is due to the fact that the array
elements used in the benchmark are 32-bit signed integers, which take up 4 bytes
each.}. For $\texttt{csize = 128} = 512$ bytes and $\texttt{stride = 1} = 4$
bytes the mean access time is $0.124$ nanoseconds, while for $\texttt{csize =
$2^{20}$} = 4$MB and for $\texttt{stride = $2^{19}$} = 2$MB the mean access time
is $1.156$ nanoseconds. The first set of parameters performs well thanks to the
low \textit{stride} value, thus achieving very good space locality and
maximizing cache hits. However, the second set of parameters achieves good
performance as well thanks to the few values accessed with each pass, thus
improving the temporal locality of each address accessed. This observation
applies for the few last data points in each data series of Figure
\ref{fig:mem}, i.e. for \textit{stride} values close to \textit{csize}. 

\subsection{Analyzing Benchmark Results}

\begin{figure}[t]
    \begin{subfigure}{0.5\textwidth}
\includegraphics[width=\textwidth]{generic_macos.pdf}
        \caption{Personal laptop}
        \label{fig:mem:laptop}
    \end{subfigure}
    \begin{subfigure}{0.5\textwidth}
\includegraphics[width=\textwidth]{generic_cluster.pdf}
        \caption{Cluster}
        \label{fig:mem:cluster}
    \end{subfigure}
    \caption{Results of the \texttt{membench.c} benchmark for both my personal
    laptop (in Figure \ref{fig:mem:laptop}) and the cluster (in Figure
    \ref{fig:mem:cluster}).}
    \label{fig:mem}
\end{figure}

The \texttt{membench.c} benchmark results for my personal laptop (Macbook Pro
2018 with a Core i7-8750H CPU) and the cluster are shown in figure
\ref{fig:mem}.

The memory access graph for the cluster's benchmark results shows that temporal
locality is best for small array sizes and for small \texttt{stride} values.
In particular, for array memory sizes of 16MB or lower (\texttt{csize} of $4
\cdot 2^{20}$ or lower) and \texttt{stride} values of 2048 or lower the mean
read+write time is less than 10 nanoseconds. Temporal locality is worst for
large sizes and strides, although the largest values of \texttt{stride} for each
size (like \texttt{csize / 2} and \texttt{csize / 4}) achieve better mean times
for the aformentioned effect of having \textit{stride} values close to
\textit{csize}.

The pattern that can be read from the graphs, especially the one for the
cluster, shows that the \textit{stride} axis is divided in regions showing
memory access time of similar magnitude. The boundary between the first and the
second region is a \textit{stride} value of rougly 2KB, while a \textit{stride}
of 512KB roughly separates the second and the third region. The difference in
performance between regions and the similarity of performance within regions
suggest the threshold stride values are related to changes in the use of the
cache hierarchy. In particular, the first region may characterize regions where
the L1 cache, the fastest non-register memory available, is predominantly used.
Then the second region might overlap with a more intense use of the L2 cache and
likewise between the third region and the L3 cache.
\marginpar[right text]{\color{white}\url{https://youtu.be/JzJlzGaQFoc}}

\section{Optimize Square Matrix-Matrix Multiplication  \punkte{60}}


\begin{figure}[t]
\begin{verbatim}
INPUT: A (n by n), B (n by n), n
OUTPUT: C (n by n)

s := 26 # block dimension

A_row := <matrix A converted in row major form>
C_temp := <empty s by s matrix>

for i := 0 to n by s:
    i_next := min(i + s, n)

    for j := 0 to n by s:
        j_next := min(j + s, n)

        <set all cells in C_temp to 0>

        for k := 0 to n by s:
            k_next := min(k + s, n)

            # Perform naive matrix multiplication, incrementing cells of C_temp
            # with each multiplication result
            naivemm(A_row[i, k][i_next, k_next], B[k, j][k_next, j_next],
                    C_temp[0, 0][i_next - i, j_next - j])
        end for

        C[i, j][i_next, j_next] = C_temp[0, 0][i_next - i, j_next - j]
    end for
end for
\end{verbatim}
    \caption{Pseudocode listing of my blocked matrix multiplication
    implementation. Matrix indices start from 0 (i.e. row $0$ and column $0$
    denotes the top-left-most cell in a matrix). \\ \texttt{M[a, b][c, d]} denotes
    a rectangular region of the matrix $M$ whose top-left-most cell is the cell
    in $M$ at row $a$ and column $b$ and whose bottom-right-most cell is the 
    cell in $M$ at row $c - 1$ and column $d - 1$.}
    \label{fig:algo}
\end{figure}

The file \texttt{matmult/dgemm-blocked.c} contains a C implementation of the
blocked matrix multiplication algorithm presented in the project. A pseudocode
listing of the implementation is provided in Figure \ref{fig:algo}.

In order to achieve a correct and fast execution, my implementation:

\begin{figure}[t]
\begin{center}
\begin{tikzpicture}
    \fill[blue!60!white] (4,0) rectangle (5,-2);
    \fill[blue!40!white] (4,-2) rectangle (5,-4);
    \fill[blue!60!white] (0,-4) rectangle (2,-5);
    \fill[blue!40!white] (2,-4) rectangle (4,-5);
    \fill[green!40!white] (4,-4) rectangle (5,-5);
    \draw[step=1,gray,very thin] (0,0) grid (5,-5);
    \draw[step=2] (0,0) grid (5,-5);
    \draw[step=5] (0,0) grid (5,-5);
\end{tikzpicture}
\end{center}
    \caption{Result of the block division process of a square matrix of size 5
    using a block size of 2. The 2-by-1 and 1-by-2 rectangular remainders are
    shown in blue and the square matrix of remainder size (i.e. 1) is shown in
    green.}
    \label{fig:matrix}
\end{figure}

\begin{itemize}
    \item Handles the edge cases related to the ``remainders'' in the matrix
        block division, i.e. when the division between the size of the matrix
        and the block size yields a remainder. Assuming only squared matrices
        are multiplied through the algorithm (as in the test suite provided) the
        block division could yield rectangular matrix blocks located in the last
        rows and columns of each matrix, and the bottom-right corner of the
        matrix will be contained in a square matrix block of the size of the
        remainder. The result of this process is shown in Figure
        \ref{fig:matrix};
    \item Converts matrix A into row major format. As shown in Figure
        \ref{fig:iter}, by having A in row major format and B in column major
        format, iterations across matrix block in the inner most loop of the
        algorithm (the one calling \textit{naivemm}) cache hits are maximised by
        achieving space locality between the blocks used. This achieved
        approximately an increase of performance of two percentage points in
        terms of CPU utilization (i.e. from a baseline of $4\%$ to $6\%$),
    \item Caches the result of each innermost iteration into a temporary matrix
        of block size before storing it into matrix C. This achieves better
        space locality when \textit{naivemm} needs to store values in matrix C.
        The block size temporary matrix has virtually no stride and thus cache
        hits are maximised. The copy operation is implemented with bulk copy
        \texttt{memcpy} calls. This optimization achieves an extra half of a
        percentage point in terms of CPU utilization (i.e. from the $6\%$
        discussed above to roughly $6.5\%$).
    \item Exploits some compiler optimizations, namely using the compiler
        optimizer at the \texttt{-O3} setting and using the \texttt{-ffast-math}
        and \texttt{-march=haswell} to respectively apply some floating point
        arithmetic optimizations and to set the compiler target to the exact ISA
        of the cluster's processor. Note that these flags are applied to all
        implementations used in the benchmark as the flags were added in the 
        \texttt{Makefile}. In addition, the \texttt{naivemm} algoritms was
        inlined in the actual \texttt{dgemm-blocked.c} source code instead of
        being separated into a function to be called to achieve better
        performance and compiler optimizations. All these changes increased the
        CPU utilization from $6.5\%$ to an average of $13.45\%$.
\end{itemize}

The chosen matrix block size for running the benchmark on the cluster is:

$$s = 32$$

as shown in the pseudocode. This value has been obtained by running an empirical
binary search on the value using the benchmark as a metric, i.e. by running
\texttt{./run\_matrixmult.sh} several times with different values. For square
blocks (i.e. the worst case) the total size for the matrix $A$ and $B$ sub-block
and the \texttt{C\_temp} temporary matrix block for $C$ is:

$$\mathrm{Bytes} = \mathrm{cellSize} * s^2 * 3 = 8 * 32^2 * 3 = 24576$$

given that a double-precision floating point number, the data type used for
matrix cells in the scope of this project, is 8 bytes long. The obtained total
bytes size is fairly close to the L1 cache size of the processor used in the
cluster ($32\mathrm{Kb} = 32768$ bytes), which is expected given that the
algorithm needs to exploit fast memory as much as possible. The reason the
empirically best value results in a theoretical cache allocation that is only
half of the complete L1 cache size is due to some real-life factors. For example,
cache misses tipically result in aligned page loads which may load unnecessary
data.

A potential way to exploit the different cache levels is to apply the blocked
matrix algorithm iteratively multiple times. For example, OpenBLAS implements
DGEMM by having two levels of matrix blocks to better exploit the L2 and L3
caches found on most processors.

\begin{figure}[t]
\begin{center}
\begin{tikzpicture}
    \node[align=center] at (2.5,0.5) {Matrix A};
    \fill[orange!10!white] (0,0) rectangle (2,-2);
    \fill[orange!25!white] (2,0) rectangle (4,-2);
    \fill[orange!40!white] (4,0) rectangle (5,-2);
    
    \draw[step=1,gray,very thin] (0,0) grid (5,-5);
    \draw[step=2,black,thick] (0,0) grid (5,-5);
    \draw[step=5,black,thick] (0,0) grid (5,-5);
    
    \draw[-to,step=1,red,very thick] (0.5,-0.5) -- (4.5,-0.5); 
    \draw[-to,step=1,red,very thick] (0.5,-1.5) -- (4.5,-1.5); 
    \draw[-to,step=1,red,very thick] (0.5,-2.5) -- (4.5,-2.5); 
    \draw[-to,step=1,red,very thick] (0.5,-3.5) -- (4.5,-3.5); 
    \draw[-to,step=1,red,very thick] (0.5,-4.5) -- (4.5,-4.5); 
    
    \node[align=center] at (8.5,0.5) {Matrix B};
    \fill[orange!10!white] (6,0) rectangle (8,-2);
    \fill[orange!25!white] (6,-2) rectangle (8,-4);
    \fill[orange!40!white] (6,-4) rectangle (8,-5);
    
    \draw[step=1,gray,very thin] (6,0) grid (11,-5);
    \draw[step=2,black,thick] (6,0) grid (11,-5);
    \draw[step=5,black,thick] (6,0) grid (11,-5);
    \draw[black,thick] (11,0) -- (11,-5);

    \draw[-to,step=1,red,very thick] (6.5,-0.5) -- (6.5,-4.5); 
    \draw[-to,step=1,red,very thick] (7.5,-0.5) -- (7.5,-4.5); 
    \draw[-to,step=1,red,very thick] (8.5,-0.5) -- (8.5,-4.5); 
    \draw[-to,step=1,red,very thick] (9.5,-0.5) -- (9.5,-4.5); 
    \draw[-to,step=1,red,very thick] (10.5,-0.5) -- (10.5,-4.5); 
\end{tikzpicture}
\end{center}
    \caption{Inner most loop iteration of the blocked GEMM algorithm across
    matrices A and B. The red lines represent the ``majorness'' of each matrix
    (A is converted by the algorithm in row-major form, while B is given and
    used in column-major form). The shades of orange represent the blocks used
    in each iteration.}
    \label{fig:iter}
\end{figure}

The results of the matrix multiplication benchmark for the naive, blocked, and
BLAS implementations are shown in Figure \ref{fig:bench} as a graph of GFlop/s
over matrix size or in Figure \ref{fig:benchtab} as a table. The blocked
implementation achieves up to 200\% more FLOPS than the naive implementation for
the largest matrix dimensions. However, the blocked implementation achives
roughly an eighth of the FLOPS the Intel MKL BLAS based implementation achieves.

I was unable to run this benchmark suite on my personal machine due to Intel MKL
installation issues that prevented the code to compile.

\begin{figure}[t]
    \includegraphics[width=\textwidth]{timing.pdf}
    \caption{GFlop/s per matrix size of the matrix multiplication benchmark for the naive,
    blocked, and BLAS implementations. The Y-axis is log-scaled.}
    \label{fig:bench}
\end{figure}

\begin{figure}[t]
\begin{center}
\begin{tabular}{c|cc|cc|cc}
    \toprule
    & \multicolumn{2}{c|}{Naive} & \multicolumn{2}{c|}{Blocked} &
    \multicolumn{2}{c}{BLAS} \\
    \makecell{Size} & \makecell{MFLOPS} &
    \makecell{\% CPU} & \makecell{MFLOPS} &
    \makecell{\% CPU} & \makecell{MFLOPS} &
    \makecell{\% CPU} \\
    \midrule
31	 & 3140.45	 & 8.53 & 3844.56	 & 10.45 & 25677.4	 & 69.78 \\
32	 & 3364.78	 & 9.14 & 5342.55	 & 14.52 & 28952.1	 & 78.67 \\
96	 & 2703.08	 & 7.35 & 5620.08	 & 15.27 & 32816.4	 & 89.18 \\
97	 & 2729.68	 & 7.42 & 4754.1	 & 12.92 & 31699.2	 & 86.14 \\
127	 & 2556.58	 & 6.95 & 4977.82	 & 13.53 & 30274.5	 & 82.27 \\
128	 & 1803.41	 & 4.90 & 4817.8	 & 13.09 & 32721.7	 & 88.92 \\
129	 & 2669.26	 & 7.25 & 4594.25	 & 12.48 & 31746.4	 & 86.27 \\
191	 & 2290.09	 & 6.22 & 4931.27	 & 13.40 & 32263.1	 & 87.67 \\
192	 & 1801.66	 & 4.90 & 5549.67	 & 15.08 & 35491.2	 & 96.44 \\
229	 & 2218.61	 & 6.03 & 4982.59	 & 13.54 & 34557.2	 & 93.91 \\
255	 & 2178.15	 & 5.92 & 4528.43	 & 12.31 & 33771.3	 & 91.77 \\
256	 & 808.413	 & 2.20 & 4652.68	 & 12.64 & 35221.1	 & 95.71 \\
257	 & 2238.93	 & 6.08 & 4512.33	 & 12.26 & 33807.9	 & 91.87 \\
319	 & 2174.45	 & 5.91 & 5093.38	 & 13.84 & 34415.8	 & 93.52 \\
320	 & 1612.13	 & 4.38 & 5674.61	 & 15.42 & 36500.2	 & 99.19 \\
321	 & 2173.64	 & 5.91 & 5111.09	 & 13.89 & 35508.1	 & 96.49 \\
417	 & 2125.36	 & 5.78 & 5143.98	 & 13.98 & 36157.6	 & 98.25 \\
479	 & 2107.13	 & 5.73 & 5152.51	 & 14.00 & 36186.4	 & 98.33 \\
480	 & 1848.43	 & 5.02 & 5703	 & 15.50 & 37971.3	 & 103.18 \\
511	 & 2112.99	 & 5.74 & 4479.96	 & 12.17 & 35144	 & 95.50 \\
512	 & 801.127	 & 2.18 & 4596.26	 & 12.49 & 37362.5	 & 101.53 \\
639	 & 1881.94	 & 5.11 & 5168.59	 & 14.05 & 36989.1	 & 100.51 \\
640	 & 815.847	 & 2.22 & 5232.97	 & 14.22 & 38267.8	 & 103.99 \\
767	 & 1825.75	 & 4.96 & 4701.09	 & 12.77 & 37220.8	 & 101.14 \\
768	 & 812.933	 & 2.21 & 4826.12	 & 13.11 & 38744	 & 105.28 \\
769	 & 1825.38	 & 4.96 & 4686.21	 & 12.73 & 37076.1	 & 100.75 \\
    \bottomrule
\end{tabular}
\end{center}
    \caption{MFlop/s and CPU utlisation per matrix size of the matrix
    multiplication benchmark for the naive, blocked, and BLAS implementations.}
    \label{fig:benchtab}
\end{figure}
\end{document}