new

diff --git a/BookGPU/BookGPU.tex b/BookGPU/BookGPU.tex
index d96b4e9db4e767a71129d1f1d200c8f5a58d16e4..e57e5958a06c50228b05ea514e57f31dc0c76380 100755 (executable)
--- a/BookGPU/BookGPU.tex
+++ b/BookGPU/BookGPU.tex
@@ -154,7 +154,7 @@
 
 
 \makeindex
 
 
 \makeindex

-%\includeonly{Chapters/chapter10/ch10}
+%\includeonly{Chapters/chapter19/ch19}

 \DeclareCaptionLabelSeparator{colon}{. }
 \begin{document}
 
 \DeclareCaptionLabelSeparator{colon}{. }
 \begin{document}
 

diff --git a/BookGPU/Chapters/chapter18/ch18.tex b/BookGPU/Chapters/chapter18/ch18.tex
index 71d2b84ee02307ade24a3af9a91df0c9e0bb82e8..d161c762743db39f988c5184b14a94ae4f4b22b0 100755 (executable)
--- a/BookGPU/Chapters/chapter18/ch18.tex
+++ b/BookGPU/Chapters/chapter18/ch18.tex
@@ -247,26 +247,7 @@ satisfies the Devaney's definition of chaos.
 
 
 
 
 
 

-\lstset{language=C,caption={C code of the sequential PRNG based on chaotic iterations},label={algo:seqCIPRNG}}
-\begin{small}
-\begin{lstlisting}
-
-unsigned int CIPRNG() {
-  static unsigned int x = 123123123;
-  unsigned long t1 = xorshift();
-  unsigned long t2 = xor128();
-  unsigned long t3 = xorwow();
-  x = x^(unsigned int)t1;
-  x = x^(unsigned int)(t2>>32);
-  x = x^(unsigned int)(t3>>32);
-  x = x^(unsigned int)t2;
-  x = x^(unsigned int)(t1>>32);
-  x = x^(unsigned int)t3;
-  return x;
-}
-\end{lstlisting}
-\end{small}
-
+\lstinputlisting[label=algo:seqCIPRNG,caption={C code of the sequential PRNG based on chaotic iterations}]{Chapters/chapter18/code2.cu}

 
 
 In Listing~\ref{algo:seqCIPRNG} a sequential  version of the proposed PRNG based
 
 
 In Listing~\ref{algo:seqCIPRNG} a sequential  version of the proposed PRNG based


@@ -298,8 +279,7 @@ simultaneously. In general,  the larger the number of  threads is, the
 more local  memory is  used, and the  less branching  instructions are
 used  (if,  while,  etc.) and so,  the  better the  performances  on  GPU  are.
 Obviously, having these requirements in  mind, it is possible to build
 more local  memory is  used, and the  less branching  instructions are
 used  (if,  while,  etc.) and so,  the  better the  performances  on  GPU  are.
 Obviously, having these requirements in  mind, it is possible to build

-a   program    similar   to    the   one   presented    in  Listing 
-\ref{algo:seqCIPRNG}, which computes  pseudorandom numbers on GPU.  To
+a   program    similar   to    the   one   presented    in  Listing~\ref{algo:seqCIPRNG}, which computes  pseudorandom numbers on GPU.  To

 do  so,  we  must   first  recall  that  in  the  CUDA~\cite{Nvid10}
 environment,    threads    have     a    local    identifier    called
 \texttt{ThreadIdx},  which   is  relative  to   the  block  containing
 do  so,  we  must   first  recall  that  in  the  CUDA~\cite{Nvid10}
 environment,    threads    have     a    local    identifier    called
 \texttt{ThreadIdx},  which   is  relative  to   the  block  containing


@@ -337,7 +317,7 @@ NumThreads: number of threads\;}
 \If{threadIdx is concerned by the computation} {
   retrieve data from InternalVarXorLikeArray[threadIdx] in local variables\;
   \For{i=1 to n} {
 \If{threadIdx is concerned by the computation} {
   retrieve data from InternalVarXorLikeArray[threadIdx] in local variables\;
   \For{i=1 to n} {

-    compute a new PRNG as in Listing\ref{algo:seqCIPRNG}\;
+    compute a new PRNG as in Listing~\ref{algo:seqCIPRNG}\;

     store the new PRNG in NewNb[NumThreads*threadIdx+i]\;
   }
   store internal variables in InternalVarXorLikeArray[threadIdx]\;
     store the new PRNG in NewNb[NumThreads*threadIdx+i]\;
   }
   store internal variables in InternalVarXorLikeArray[threadIdx]\;

diff --git a/BookGPU/Chapters/chapter18/code2.cu b/BookGPU/Chapters/chapter18/code2.cu
new file mode 100644 (file)
index 0000000..0eb74bd
--- /dev/null
+++ b/BookGPU/Chapters/chapter18/code2.cu
@@ -0,0 +1,14 @@
+{
+unsigned int CIPRNG() {
+  static unsigned int x = 123123123;
+  unsigned long t1 = xorshift();
+  unsigned long t2 = xor128();
+  unsigned long t3 = xorwow();
+  x = x^(unsigned int)t1;
+  x = x^(unsigned int)(t2>>32);
+  x = x^(unsigned int)(t3>>32);
+  x = x^(unsigned int)t2;
+  x = x^(unsigned int)(t1>>32);
+  x = x^(unsigned int)t3;
+  return x;
+}

diff --git a/BookGPU/Chapters/chapter19/ch19.tex b/BookGPU/Chapters/chapter19/ch19.tex
index fdf533d46bb0d01853770689c12ec0ccea95e697..480e5287722dc2512f362a7e45e3dde375fc1f49 100755 (executable)
--- a/BookGPU/Chapters/chapter19/ch19.tex
+++ b/BookGPU/Chapters/chapter19/ch19.tex
@@ -7,50 +7,50 @@
 \section{Introduction}
 \label{ch19:intro}
 
 \section{Introduction}
 \label{ch19:intro}
 

-The Number Field Sieve (NFS)\index{iterative methods!Number Field Sieve} is the current state-of-the-art integer factorization method. It requires the solution of a large sparse linear system over GF(2) (called the linear algebra step). The Block Wiedemann\index{Number Field Sieve!Block Wiedemann} (BW)\cite{ch19:bw} algorithm can be used to solve such a large sparse linear system efficiently using iterative sparse matrix vector multiplication (SpMV).
+The Number Field Sieve (NFS)\index{iterative methods!Number Field Sieve} is the current state-of-the-art integer factorization method. It requires the solution of a large sparse linear system over Galois Field GF(2) (called the linear algebra step). The Block Wiedemann\index{Number Field Sieve!Block Wiedemann} (BW)\cite{ch19:bw} algorithm can be used to solve such a large sparse linear system efficiently using iterative sparse matrix vector multiplication (SpMV).

 
 Recent integer factorization efforts have been using CPU clusters to solve the large sparse linear system \cite{ch19:kilobit,ch19:rsa768}. The RSA-768 factorization \cite{ch19:rsa768}, for example, reported a runtime of 3 months for the linear algebra step on a cluster with 48 AMD dual hex-core CPUs. Previous work on parallelizing the linear algebra step focused on using CPU clusters and grids \cite{ch19:aoki,ch19:hwang,ch19:grid,ch19:hetero768}. In this chapter, we present a CUDA approach that can be used to accelerate the costly iterative SpMV operation for matrices derived from NFS.
 
 
 Recent integer factorization efforts have been using CPU clusters to solve the large sparse linear system \cite{ch19:kilobit,ch19:rsa768}. The RSA-768 factorization \cite{ch19:rsa768}, for example, reported a runtime of 3 months for the linear algebra step on a cluster with 48 AMD dual hex-core CPUs. Previous work on parallelizing the linear algebra step focused on using CPU clusters and grids \cite{ch19:aoki,ch19:hwang,ch19:grid,ch19:hetero768}. In this chapter, we present a CUDA approach that can be used to accelerate the costly iterative SpMV operation for matrices derived from NFS.
 

-The memory access pattern in the SpMV operation generally consists of regular access patterns over the matrix and irregular access patterns over the vector. The irregular access pattern over the vector is a challenge that is pronounced more on the GPU than on the CPU, because of the smaller cache and the restrictive memory access pattern requirement to achieve maximum performance. However, a high-end GPU has an order-of-magnitude higher bandwidth than a high-end CPU; e.g. a GeForce GTX 580 has 192.4 GB/s memory bandwidth, while an Intel Core-i7 has a maximum of 25.6 GB/s memory bandwidth. 
+The memory access pattern in the SpMV operation generally consists of regular access patterns over the matrix and irregular access patterns over the vector. The irregular access pattern over the vector is a challenge that is more pronounced  on the GPU than on the CPU, because of the smaller cache and the restrictive memory access pattern requirement to achieve maximum performance. However, a high-end GPU has an order-of-magnitude higher bandwidth than a high-end CPU; e.g., a GeForce GTX 580 has 192.4 GB/s memory bandwidth, while an Intel Core-i7 has a maximum of 25.6 GB/s memory bandwidth. 

 
 

-SpMV on the GPU has been explored previously in several papers \cite{ch19:nvidia,ch19:exactspmv,ch19:bellpack,ch19:sle} for matrices derived from scientific computing applications. However, sparse matrices derived from NFS have generally different properties, i.e. they are larger, have a few dense rows and have many extremely sparse rows. The large size of the matrix causes the BL and BW algorithms to require a large number of SpMV iterations. This means that the time spent for matrix preprocessing and matrix data transfer to the GPU memory are negligible compared to the total runtime. Thus, approaches to SpMV on GPUs for NFS matrices may be different from previously published GPU SpMV approaches. We further present an extension of our SpMV method for binary-valued NFS matrices to single-precision floating-point matrices.
+SpMV on the GPU has been explored previously in several papers \cite{ch19:nvidia,ch19:exactspmv,ch19:bellpack,ch19:sle} for matrices derived from scientific computing applications. However, sparse matrices derived from NFS have generally different properties, i.e. they are larger, have a few dense rows, and have many extremely sparse rows. The large size of the matrix causes the BL and BW algorithms to require a large number of SpMV iterations. This means that the time spent for matrix preprocessing and matrix data transfer to the GPU memory is negligible compared to the total runtime. Thus, approaches to SpMV on GPUs for NFS matrices may be different from previously published GPU SpMV approaches. We further present an extension of our SpMV method for binary-valued NFS matrices to single-precision floating-point matrices.

 
 

-\section{Block Wiedemann Algorithm}
+\section{Block Wiedemann algorithm}

 \label{ch19:block-wiedemann}
 
 \label{ch19:block-wiedemann}
 

-The BW algorithm heuristically finds $n$ vectors in the kernel space \index{Number Field Sieve!kernel space}of a $d \times d$ binary matrix $B$, $n$ is one of two parameters $m, n$, called blocking factors\index{Number Field Sieve!blocking factors}. BW consists of the following steps:
+The BW algorithm heuristically finds $n$ vectors in the kernel space \index{Number Field Sieve!kernel space}of a $d \times d$ binary matrix $B$; $n$ is one of two parameters $m, n$, called blocking factors\index{Number Field Sieve!blocking factors}. BW consists of the following steps:

 
 \begin{itemize}
 \item \textbf{Step 1 (BW1):} Compute the matrix sequence 
 
 \begin{equation}
 
 \begin{itemize}
 \item \textbf{Step 1 (BW1):} Compute the matrix sequence 
 
 \begin{equation}

-A_i = x \cdot B^i \cdot y, \forall i=1,...,{\frac{d}{m}}+{ \frac{d}{n} }+O(1)
+A_i = x \cdot B^i \cdot y, \forall i=1,...,{\frac{d}{m}}+{ \frac{d}{n} }+O(1),

 \end{equation}
 
 \end{equation}
 

- where $x,y$ are randomly chosen binary matrices of size $m \times d$ and $d \times n$ respectively.
+ where $x,y$ are randomly chosen binary matrices of size $m \times d$ and $d \times n$, respectively.

 \item \textbf{Step 2 (BW2):} The Berlekamp-Massey\index{Number Field Sieve!Berlekamp-Massey} algorithm \cite{ch19:Thome:subqad} is used to compute a generating polynomial of the matrix sequence $A$ from BW1 in the form
 \begin{equation}
 \item \textbf{Step 2 (BW2):} The Berlekamp-Massey\index{Number Field Sieve!Berlekamp-Massey} algorithm \cite{ch19:Thome:subqad} is used to compute a generating polynomial of the matrix sequence $A$ from BW1 in the form
 \begin{equation}

-F(X)= \sum_{i=1}^{{ \frac{d}{n} }+O(1)} {C_i \cdot X^i}
+F(X)= \sum_{i=1}^{{ \frac{d}{n} }+O(1)} {C_i \cdot X^i},

 \end{equation}
 \end{equation}

-where $C_i$ is a $n \times n$ binary matrix
+where $C_i$ is an $n \times n$ binary matrix

 
 \item \textbf{Step 3 (BW3):} Compute the sequence of matrices $S$ of size $d \times n$ such that
 \begin{equation}
 
 \item \textbf{Step 3 (BW3):} Compute the sequence of matrices $S$ of size $d \times n$ such that
 \begin{equation}

-S_i = B^i \cdot y \cdot C_i, \forall i=1,...,{\frac{d}{n}}+O(1)
+S_i = B^i \cdot y \cdot C_i, \forall i=1,...,{\frac{d}{n}}+O(1).

 \end{equation}
 
 \end{equation}
 

-With high probability, $B \cdot \sum{S_i} = 0$ for which $\sum{S_i}$ is output as a solution.
+With high probability, $B \cdot \sum{S_i} = 0$, for which $\sum{S_i}$ is output as a solution.

 \end{itemize}
 
 \end{itemize}
 

-We can treat $x$, $y$ and $S_i$ as vectors of block width $m$ or $n$. Assuming $B$ is a sparse matrix with $\gamma$ non-zeros per row on average and ignoring $m$ and $n$ as constant, the complexity of BW1 and BW3 is $O(\gamma d^2)$. Using the subquadratic algorithm by Thome \cite{ch19:Thome:subqad}, BW2 has a complexity of $O(d \log^2{d})$. Thus our approach to accelerate BW is based on accelerating BW1 and BW3 on a GPU while BW2 is still done on a CPU.
+We can treat $x$, $y$, and $S_i$ as vectors of block width $m$ or $n$. Assuming $B$ is a sparse matrix with $\gamma$ nonzeros per row on average and ignoring $m$ and $n$ as constant, the complexity of BW1 and BW3 is $O(\gamma d^2)$. Using the subquadratic algorithm by Thomé \cite{ch19:Thome:subqad}, BW2 has a complexity of $O(d \log^2{d})$. Thus our approach to accelerate BW is based on accelerating BW1 and BW3 on a GPU while BW2 is still done on a CPU.

 
 

-\section{SpMV OVER GF(2) for NFS matrices using existing formats on GPUs}
+\section{SpMV over GF(2) for NFS matrices using existing formats on GPUs}

 \label{ch19:existingspmv}
        In this section, we review a few relevant previously published sparse matrix formats on GPUs and study their performance when applied to sparse matrices over GF(2) derived from integer factorization with NFS.
        
 \label{ch19:existingspmv}
        In this section, we review a few relevant previously published sparse matrix formats on GPUs and study their performance when applied to sparse matrices over GF(2) derived from integer factorization with NFS.
        

-       We consider the sparse binary matrix $B$ of size $d \times d$ and a dense vector $X$ of size $d \times n$ bit, where $n$ is called the \emph{blocking factor}. Mordern processors generally support 64-bit operations. Thus, typical blocking factors are of the form of $64 \cdot k, \; k \in \mathbb{N}$. Note that doubling the blocking factor roughly halves the number of SpMV iterations required but doubles the input vector size. 
+       In Algorithm~\ref{ch19:algo1}, we consider the sparse binary matrix $B$ of size $d \times d$ and a dense vector $X$ of size $d \times n$ bit, where $n$ is called the \emph{blocking factor}. Modern processors generally support 64-bit operations. Thus, typical blocking factors are of the form of $64 \cdot k, \; k \in \mathbb{N}$. Note that doubling the blocking factor roughly halves the number of SpMV iterations required but doubles the input vector size. 

 
 

-       For all $0 \leq i \leq d-1$ let $c\_index[i]$ column index array of $B[i]$ which contains the indices of the non-zero entries of row $i$. Then, the following pseudocode shows a single SpMV iteration of $B$ with input vector $X$ and result vector $Y$.
+       For all $0 \leq i \leq d-1$ let $c\_index[i]$ column index array of $B[i]$ which contains the indices of the nonzero entries of row $i$. Then, the following pseudocode shows a single SpMV iteration of $B$ with input vector $X$ and result vector $Y$.

 
 %% \begin{algorithm}
 %% \caption{SpMV-COLUMN-INDEX}
 
 %% \begin{algorithm}
 %% \caption{SpMV-COLUMN-INDEX}


@@ -70,7 +70,8 @@ We can treat $x$, $y$ and $S_i$ as vectors of block width $m$ or $n$. Assuming $
 
 
 \begin{algorithm}
 
 
 \begin{algorithm}

-\caption{SpMV-COLUMN-INDEX}
+\caption{SpMV-column-index}
+\label{ch19:algo1}

 %      \begin{algorithmic}
        \KwIn{$c\_index[i]$: column array of rows $i$ of $B$}
        \KwIn{ $X$: input vector}
 %      \begin{algorithmic}
        \KwIn{$c\_index[i]$: column array of rows $i$ of $B$}
        \KwIn{ $X$: input vector}


@@ -79,7 +80,7 @@ We can treat $x$, $y$ and $S_i$ as vectors of block width $m$ or $n$. Assuming $
        $Y[i]=0$\;   
    \ForAll {$ind$ $\in$ $c\_index[i]$} {
                $Y[i] = Y[i] \oplus X[ind]$\;
        $Y[i]=0$\;   
    \ForAll {$ind$ $\in$ $c\_index[i]$} {
                $Y[i] = Y[i] \oplus X[ind]$\;

-               \CommentSty{$\oplus:$ bitwise XOR operation}
+               \CommentSty{$\oplus:$ bitwise xor operation}

 }
 }
 %      \end{algorithmic}
 }
 }
 %      \end{algorithmic}


@@ -88,19 +89,19 @@ We can treat $x$, $y$ and $S_i$ as vectors of block width $m$ or $n$. Assuming $
        
        The costly operations in the SpMV pseudocode are the memory accesses for loading $ind$, $X[ind]$, $Y[i]$ and storing $Y[i]$. To speed up those operations on any architecture a common approach is to design a cache-friendly order of accessing the memory. The order is especially important for the vectors $X$ and $Y$, since their memory locations might be accessed multiple times. Accesses to $Y$ can be minimized by storing the intermediate results in fast memory and only write the accumulated result to $Y$.
        
        
        The costly operations in the SpMV pseudocode are the memory accesses for loading $ind$, $X[ind]$, $Y[i]$ and storing $Y[i]$. To speed up those operations on any architecture a common approach is to design a cache-friendly order of accessing the memory. The order is especially important for the vectors $X$ and $Y$, since their memory locations might be accessed multiple times. Accesses to $Y$ can be minimized by storing the intermediate results in fast memory and only write the accumulated result to $Y$.
        

-       CUDA implementations of SpMV generally store both matrix and vectors in \emph{global memory}. Memory accesses to $B$ and $Y$ can usually be coalesced. Memory accesses to $X$ are random and non-coalesced thus having higher latency, but \emph{texture memory} can be used to take advantage of the \emph{texture cache}. Intermediate results could utilize \emph{shared memory} which has low latency however is very limited in size and only shared between threads of the same block. Shared memory is further divided into \emph{banks}. When threads in a warp access the shared memory, \emph{bank conflicts} should be avoided otherwise the access will be serialized.
+       CUDA implementations of SpMV generally store both matrix and vectors in \emph{global memory}. Memory accesses to $B$ and $Y$ can usually be coalesced. Memory accesses to $X$ are random and noncoalesced, thus having higher latency, but \emph{texture memory} can be used to take advantage of the \emph{texture cache}. Intermediate results could utilize \emph{shared memory}, which has low latency but is very limited in size and only shared between threads of the same block. Shared memory is further divided into \emph{banks}. When threads in a warp access the shared memory, \emph{bank conflicts} should be avoided, otherwise the access will be serialized.

        
        

-       We now briefly review the CUDA implementation of a number of SpMV matrix formats published in previous papers \cite{ch19:nvidia,ch19:sle}. We also describe the data structure in which the matrix in Figure \ref{fig:ex_matrix} is stored using the format as an example. Please note that when storing a binary matrix, the $value$ array can be removed.
+       We now briefly review the CUDA implementation of a number of SpMV matrix formats published in previous papers \cite{ch19:nvidia,ch19:sle}. We also describe the data structure in which the matrix in Figure \ref{fig:ex_matrix} is stored using different formats as examples. Please note that when storing a binary matrix, the $value$ array can be removed.

        
 \begin{figure}[t]
        \centering
        \includegraphics[width=4cm]{Chapters/chapter19/fig/ex_matrix.pdf} 
        
 \begin{figure}[t]
        \centering
        \includegraphics[width=4cm]{Chapters/chapter19/fig/ex_matrix.pdf} 

-       \caption{An example square matrix of size $6 \times 6$ (zero values are not shown)}
+       \caption{An example square matrix of size $6 \times 6$ (zero values are not shown).}

        \label{fig:ex_matrix}
 \end{figure}
 
        \subsubsection*{Coordinate list (COO)\index{Compressed storage format!COO}}
        \label{fig:ex_matrix}
 \end{figure}
 
        \subsubsection*{Coordinate list (COO)\index{Compressed storage format!COO}}

-For each non-zero, both its column and row indices are explicitly stored. The Cusp implementation \cite{ch19:cusp} stores elements in sorted order of row indices ensuring that entries with the same row index are stored contiguously. 
+For each nonzero, both its column and row indices are explicitly stored. The Cusp implementation \cite{ch19:cusp} stores elements in sorted order of row indices ensuring that entries with the same row index are stored contiguously. 

 
        \begin{lstlisting}[caption={}]
        coo.row_index = {0, 1, 1, 2, 2, 3, 3,  3, 4, 5,  5}
 
        \begin{lstlisting}[caption={}]
        coo.row_index = {0, 1, 1, 2, 2, 3, 3,  3, 4, 5,  5}


@@ -108,7 +109,7 @@ For each non-zero, both its column and row indices are explicitly stored. The Cu
        coo.value     = {3, 1, 5, 2, 4, 6, 8, 10, 9, 7, 11}
        \end{lstlisting}
 
        coo.value     = {3, 1, 5, 2, 4, 6, 8, 10, 9, 7, 11}
        \end{lstlisting}
 

-       \subsubsection*{Compressed Sparse Row (CSR)\index{Compressed storage format!CSR}} Non-zeros are sorted by the row index, and only their column indices are explicitly stored in a column array. Additionally, the vector $row\_start$ stores indices of the first non-zero element of each row in the column array.
+       \subsubsection*{Compressed sparse row (CSR)\index{Compressed storage format!CSR}} Nonzeros are sorted by the row index, and only their column indices are explicitly stored in a column array. Additionally, the vector $row\_start$ stores indices of the first nonzero element of each row in the column array.

        
        \begin{lstlisting}[caption={}]
        csr.row_start = {0, 1, 3, 5, 8, 9, 12}
        
        \begin{lstlisting}[caption={}]
        csr.row_start = {0, 1, 3, 5, 8, 9, 12}


@@ -116,7 +117,7 @@ For each non-zero, both its column and row indices are explicitly stored. The Cu
        csr.value     = {3, 1, 5, 2, 4, 6, 8, 10, 9, 7, 11}
        \end{lstlisting}
 
        csr.value     = {3, 1, 5, 2, 4, 6, 8, 10, 9, 7, 11}
        \end{lstlisting}
 

-       \subsubsection*{Ellpack (ELL)\index{Compressed storage format!ELL}} Let $K$ be the maximum number of non-zero elements in any row of the matrix. Then, for each row, ELL stores exactly $K$ elements (extra padding is required for rows that contain less than $K$ non-zero elements). Only column indices are required to store in an array, the row index can be implied since exactly $K$ elements are stored per row. The Cusp implementation store the column indices in a transposed manner so that consecutive threads can access consecutive memory addresses.
+       \subsubsection*{Ellpack (ELL)\index{Compressed storage format!ELL}} Let $K$ be the maximum number of nonzero elements in any row of the matrix. Then, for each row, ELL stores exactly $K$ elements (extra padding is required for rows that contain fewer than $K$ nonzero elements). Only column indices are required to store in an array, the row index can be implied since exactly $K$ elements are stored per row. The Cusp implementation stores the column indices in a transposed manner so that consecutive threads can access consecutive memory addresses.

        
        \begin{lstlisting}[caption={}]
        ell.col_index = {
        
        \begin{lstlisting}[caption={}]
        ell.col_index = {


@@ -129,7 +130,7 @@ For each non-zero, both its column and row indices are explicitly stored. The Cu
                                        *, *, *, 10, *,  *}
        \end{lstlisting}
 
                                        *, *, *, 10, *,  *}
        \end{lstlisting}
 

-       \subsubsection*{Hybrid (HYB)\index{Compressed storage format!HYB}} The HYB format heuristically computes an value $K$ and store $K$ non-zeros per rows in the ELL format. When a row has more than $K$ non-zeros, the trailing non-zeros are stored in COO. This design decreases the storage overhead due to ELL padding elements and thus improve the overall performance.
+       \subsubsection*{Hybrid (HYB)\index{Compressed storage format!HYB}} The HYB format heuristically computes a value $K$ and stores $K$ nonzeros per rows in the ELL format. When a row has more than $K$ non-zeros, the trailing nonzeros are stored in COO. This design decreases the storage overhead due to ELL padding elements and thus improves the overall performance.

        \begin{lstlisting}[caption={}]
        hyb.nnz_per_row   = 2
        hyb.ell.col_index = {2, 1, 1, 0, 2, 0, *, 4, 3, 2, *,  5}
        \begin{lstlisting}[caption={}]
        hyb.nnz_per_row   = 2
        hyb.ell.col_index = {2, 1, 1, 0, 2, 0, *, 4, 3, 2, *,  5}


@@ -139,7 +140,7 @@ For each non-zero, both its column and row indices are explicitly stored. The Cu
        hyb.coo.value     = {10}
        \end{lstlisting}
 
        hyb.coo.value     = {10}
        \end{lstlisting}
 

-       \subsubsection*{Sliced Ellpack (SLE)\index{Compressed storage format!SLE}} This format partitions the matrix into horizontal slices of $S$ adjacent rows \cite{ch19:sle}. Each slice is stored in ELLPACK format. The maximum number of non-zeros may be different for each slice. An additional array $slice\_start$ is used to index the first element in each slice. The matrix rows are usually sorted by the number of non-zeros per row in order to move rows with similar number of non-zeros together. 
+       \subsubsection*{Sliced Ellpack (SLE)\index{Compressed storage format!SLE}} This format partitions the matrix into horizontal slices of $S$ adjacent rows \cite{ch19:sle}. Each slice is stored in ELLPACK format. The maximum number of nonzeros may be different for each slice. An additional array $slice\_start$ is used to index the first element in each slice. The matrix rows are usually sorted by the number of nonzeros per row in order to move rows with similar number of nonzeros together. 

        \begin{lstlisting}[caption={}]
        sle.slice_size  = 2
        sle.col_index   = {
        \begin{lstlisting}[caption={}]
        sle.slice_size  = 2
        sle.col_index   = {


@@ -153,8 +154,7 @@ For each non-zero, both its column and row indices are explicitly stored. The Cu
        sle.slice_start = {0, 4, 10, 14}
        \end{lstlisting}        
        
        sle.slice_start = {0, 4, 10, 14}
        \end{lstlisting}        
        

-       \begin{table}
-       \caption{Properties of some NFS matrices}
+       \begin{table}[htbp]

        \centering
        \begin{tabular}
        {| l || c | c | c | c | c | c |}        
        \centering
        \begin{tabular}
        {| l || c | c | c | c | c | c |}        


@@ -162,13 +162,13 @@ For each non-zero, both its column and row indices are explicitly stored. The Cu
        &   RSA-170 & RSA-190 & KILOBIT & RSA-768 \\
        \hline            
          Max dimension & 10.4M & 26.1M & 66.7M & 192M \\
        &   RSA-170 & RSA-190 & KILOBIT & RSA-768 \\
        \hline            
          Max dimension & 10.4M & 26.1M & 66.7M & 192M \\

-         Non-zeros & 994.7M & 2,451M & 9,538M & 27,797M\\
+         Nonzeros & 994.7M & 2,451M & 9,538M & 27,797M\\

          Max row weight & 5.5M & 14M & 28.2M & 82.6M\\
          Min row weight & 3 & 3 & 2 & 2\\
          Avg row weight & 95.08 & 93.6 & 143 & 144\\
        \hline
        \end{tabular}
          Max row weight & 5.5M & 14M & 28.2M & 82.6M\\
          Min row weight & 3 & 3 & 2 & 2\\
          Avg row weight & 95.08 & 93.6 & 143 & 144\\
        \hline
        \end{tabular}

-               
+       \caption{Properties of some NFS matrices}

        \label{table:rsa_matrix}        
        \end{table}
        
        \label{table:rsa_matrix}        
        \end{table}
        


@@ -177,28 +177,28 @@ The existing formats do not achieve good performance due to the special structur
 \section{A hybrid format for SpMV on GPUs}
 \label{Implementation}
 
 \section{A hybrid format for SpMV on GPUs}
 \label{Implementation}
 

-As a preprocessing step, we reorder the rows of the matrix by their \emph{row weight}, in non-increasing order. The row weight of row $j$ of $B$ is defined as the total number of non-zero elements in row $j$. We then partition the sorted matrix rows into at most four consecutive parts. Each part uses a different format. The different formats are optimized for the sparseness properties of each partition as shown in Figure \ref{fig:partitioning}. For the densest part, we use a dense format. When the matrix gets less dense, we switch to another format which we call \index{Compressed storage format!Sliced COO} Sliced COO (SCOO). SCOO has three variants, small, medium, and large. Our formats are now described in more detail.
+As a preprocessing step, we reorder the rows of the matrix by their \emph{row weight}, in nonincreasing order. The row weight of row $j$ of $B$ is defined as the total number of nonzero elements in row $j$. We then partition the sorted matrix rows into at most four consecutive parts. Each part uses a different format. The different formats are optimized for the sparseness properties of each partition as shown in Figure \ref{fig:partitioning}. For the densest part, we use a dense format. When the matrix gets less dense, we switch to another format which we call \index{Compressed storage format!Sliced COO} Sliced COO (SCOO). SCOO has three variants, small, medium, and large. Our formats are now described in more detail.

 
 \begin{figure}[t]
        \centering
        \includegraphics[scale=1.3]{Chapters/chapter19/fig/partitioning.pdf} 
 
 \begin{figure}[t]
        \centering
        \includegraphics[scale=1.3]{Chapters/chapter19/fig/partitioning.pdf} 

-       \caption{Partitioning of a row-sorted NFS matrix into four formats}
+       \caption{Partitioning of a row-sorted NFS matrix into four formats.}

        \label{fig:partitioning}
 \end{figure}
 
        \label{fig:partitioning}
 \end{figure}
 

-    \subsection{Dense Format}
-    The dense format is used for the dense part of the matrix. This format uses 1 bit per matrix entry. Within a column, 32 matrix entries are stored as a 32 bit integer. Thus, 32 rows are stored as $N$ consecutive integers.
+    \subsection{Dense format}
+    The dense format is used for the dense part of the matrix. This format uses 1 bit per matrix entry. Within a column, 32 matrix entries are stored as a 32-bit integer. Thus, 32 rows are stored as $N$ consecutive integers.

        
        

-    Each CUDA thread works on a column. Each thread fetches one element from the input vector in coalesced fashion. Then, each thread checks the 32 matrix entries one by one. When the matrix entry is a non-zero, the thread performs a XOR operation between the element from the input vector and the partial result for the row. This means each thread only accesses the input vector once to do work on up to 32 non-zeros. The partial result from each thread needs to be stored and combined to get the final result for the 32 rows. These operations are performed in CUDA shared memory.
+    Each CUDA thread works on a column. Each thread fetches one element from the input vector in coalesced fashion. Then, each thread checks the 32 matrix entries one by one. When the matrix entry is a nonzero, the thread performs an xor operation between the element from the input vector and the partial result for the row. This means each thread accesses the input vector only once to do work on up to 32 nonzeros. The partial result from each thread needs to be stored and combined to get the final result for the 32 rows. These operations are performed in CUDA shared memory.

 
 

-    The 32 threads in a warp share 32 shared memory entries to store the partial results from the 32 rows. Since all threads in a warp execute a common instruction at a time, access to these 32 entries can be made exclusive. The result from each warp in a thread block is combined using p-reduction on shared memory. The result from each thread block is combined using an atomic XOR operation on global memory.
+    The 32 threads in a warp share 32 shared memory entries to store the partial results from the 32 rows. Since all threads in a warp execute a common instruction at the same time, access to these 32 entries can be made exclusive. The result from each warp in a thread block is combined using p-reduction on shared memory. The result from each thread block is combined using an atomic xor operation on global memory.

 
 

-    When the blocking factor is larger than 64, access to the shared memory needs to be reorganized to avoid bank conflicts. Each thread can read/write up to 64 bit data at a time to the shared memory. If a thread is accessing 128 bit data for example, two read/write operations need to be performed. Thus, there will be bank conflicts if we store 128 bit data on contiguous addresses. We can avoid bank conflicts by having the threads in a warp first access consecutive 64 bit elements representing the first halves of the 128 bit elements. Then, the threads again access consecutive 64 bit elements representing the second halves of the 128 bit elements. The same modification can be applied to other formats as well.
+    When the blocking factor is larger than 64, access to the shared memory needs to be reorganized to avoid bank conflicts. Each thread can read/write up to 64-bit data at a time to the shared memory. If a thread is accessing 128-bit data for example, two read/write operations need to be performed. Thus, there will be bank conflicts if we store 128-bit data on contiguous addresses. We can avoid bank conflicts by having the threads in a warp first access consecutive 64-bit elements representing the first halves of the 128-bit elements. Then, the threads again access consecutive 64-bit elements representing the second halves of the 128-bit elements. The same modification can be applied to other formats as well.

 
        \subsection{Sliced COO}
 
        \subsection{Sliced COO}

-       The SCOO format is adapted from the CADO-NFS software for CPUs \cite{ch19:cadonfs}. The aim is to reduce irregular accesses to the input vector and increase the texture cache hit rate. SCOO stores the column index and the row index of each non-zero. A number of consecutive rows form a slice. Non-zeros within a slice are sorted by their column index. We give an example of SCOO in Figure \ref{fig:scoo-1}, threads working on a slice can access contiguous elements from the input vector. In the figure, S1, S2 are two slices. Entries in the input vector denote the corresponding non-zero elements that the input vector element is multiplied with.
+       The SCOO format is adapted from the CADO-NFS software for CPUs \cite{ch19:cadonfs}. The aim is to reduce irregular accesses to the input vector and increase the texture cache hit rate. SCOO stores the column index and the row index of each nonzero. A number of consecutive rows form a slice. Nonzeros within a slice are sorted by their column index. We give an example of SCOO in Figure \ref{fig:scoo-1}, threads working on a slice can access contiguous elements from the input vector. In the figure, S1, S2 are two slices. Entries in the input vector denote the corresponding nonzero elements with which the input vector element is multiplied.

        
        

-       For each non-zero, two bytes are used to store the column index. However, two bytes are not enough for large RSA matrices. Thus, we further divide a slice into groups. Group $i$ contains non-zeros with column index between $i \times 2^{16}$ and $(i+1) \times 2^{16} - 1$. An additional array stores the starting position of each group in the slice. 
+       For each nonzero, two bytes are used to store the column index. However, two bytes is not enough for large RSA matrices. Thus, we further divide a slice into groups. Group $i$ contains nonzeros with column index between $i \times 2^{16}$ and $(i+1) \times 2^{16} - 1$. An additional array stores the starting position of each group in the slice. 

        
        \begin{figure}[t]
        \centering
        
        \begin{figure}[t]
        \centering


@@ -210,67 +210,66 @@ As a preprocessing step, we reorder the rows of the matrix by their \emph{row we
 
 \lstinputlisting[label=ch19:lst:spmv,caption=SpMV with SCOO format]{Chapters/chapter19/code.cu}
 
 
 \lstinputlisting[label=ch19:lst:spmv,caption=SpMV with SCOO format]{Chapters/chapter19/code.cu}
 

-       One thread block works on a slice, one group at a time. Neighboring threads work on neighboring non-zeros in the group in parallel. Each thread works on more than one row of the matrix. Thus, each thread needs some storage to store the partial result and combine them with the result from the other threads to generate the final output. Each entry costs equally as the element of the vector i.e. the blocking factor in byte. Shared memory is used as intermediate storage as they have low latency however it is limited to 48 KB per SM in Fermi. The global memory is only accessed once at the end to write the final output. The CUDA kernel for SpMV with Sliced COO format is shown in Listing \ref{ch19:lst:spmv}. 
+       One thread block works on a slice, one group at a time. Neighboring threads work on neighboring nonzeros in the group in parallel. Each thread works on more than one row of the matrix. Thus, each thread needs some storage to store the partial results and combine them with the results from the other threads to generate the final output. Each entry costs equally as the element of the vector, i.e., the blocking factor in bytes. Shared memory is used as intermediate storage as it has low latency however it is limited to 48 KB per SM in Fermi. The global memory is accessed only once at the end to write the final output. The CUDA kernel for SpMV with Sliced COO format is shown in Listing \ref{ch19:lst:spmv}. 

        
        

-       Since neighboring non-zeros may or may not come from the same row (recall that we sorted them by column index), many threads may access the same entry in the shared memory if it is used to store result of the same rows. Thus, shared memory entries are either assigned exclusively to a single thread or shared by using atomic XOR operations.    Based on the way we allocate the shared memory, we further divide the Sliced COO format into three different subformats: small, medium, and large. 
+       Since neighboring nonzeros may or may not come from the same row (recall that we sorted them by column index), many threads may access the same entry in the shared memory if it is used to store result of the same rows. Thus, shared memory entries are either assigned exclusively to a single thread or shared by using atomic xor operations.     Based on the way we allocate the shared memory, we further divide the Sliced COO format into three different subformats: small, medium, and large. 

        
        \begin{table}
        
        \begin{table}

-       \caption{Sliced COO subformat comparison (\# rows per slices is based on $n=64$) }

        \centering
        \begin{tabular}
        { | c || c | c | c |}   
        \hline  
        Sliced COO subformats & Small & Medium & Large \\
        \hline            
        \centering
        \begin{tabular}
        { | c || c | c | c |}   
        \hline  
        Sliced COO subformats & Small & Medium & Large \\
        \hline            

-         Memory Sharing & No sharing & Among Warp & Among Block\\
-         Access method & Direct & Atomic XOR & Atomic XOR\\
+         Memory sharing & No sharing & Among warp & Among block\\
+         Access method & Direct & Atomic xor & Atomic xor\\

          Bank conflict & No & No & Yes\\
          \# Rows per Slice & 12 & 192 & 6144 \\
        \hline
        \end{tabular}
          Bank conflict & No & No & Yes\\
          \# Rows per Slice & 12 & 192 & 6144 \\
        \hline
        \end{tabular}

+       \caption{Sliced COO subformat comparison (\# rows per slices is based on $n=64$).}

        \label{table:scoo}
        \end{table}
 
        \label{table:scoo}
        \end{table}
 

-       \subsubsection*{Small Sliced (SS) COO}
+       \subsubsection*{Small sliced (SS) COO}

        In this subformat, each thread has one exclusive entry in shared memory to store the partial result for each row. The assignment of the shared memory is organized such that each thread in a warp accesses only one bank and there is no bank accessed by more than one thread. Thus, there is no bank conflict. A p-reduction operation on shared memory is required to combine partial results from each thread.
 
        In this subformat, each thread has one exclusive entry in shared memory to store the partial result for each row. The assignment of the shared memory is organized such that each thread in a warp accesses only one bank and there is no bank accessed by more than one thread. Thus, there is no bank conflict. A p-reduction operation on shared memory is required to combine partial results from each thread.
 

-       The maximum number of rows per slice is calculated as \emph{size of shared memory per SM in bits / (number of threads per block * blocking factor)}. We use 512 threads per block for 64 bit blocking factor which gives 12 rows, and 256 threads per block for 128 and 256 bit blocking factor, which gives 12 and 6 rows, respectively. Hence, one byte per row index is sufficient for this subformat.
+       The maximum number of rows per slice is calculated as \emph{size of shared memory per SM in bits} / \emph{(number of threads per block * blocking factor)}. We use 512 threads per block for 64-bit blocking factor which gives 12 rows, and 256 threads per block for 128 and 256-bit blocking factor, which gives 12 and 6 rows, respectively. Hence, one byte per row index is sufficient for this subformat.

        
        

-       \subsubsection*{Medium Sliced (MS) COO}
-       In this subformat, each thread in a warp gets an entry in the shared memory to store the partial result for each row. However, this entry is shared with the threads in other warps. Access to the shared memory uses an atomic XOR operation. Each thread in a warp accesses only one bank, avoiding bank conflicts. A p-reduction operation on shared memory is required to combine the 32 partial results.
+       \subsubsection*{Medium sliced (MS) COO}
+       In this subformat, each thread in a warp gets an entry in the shared memory to store the partial result for each row. However, this entry is shared with the threads in other warps. Access to the shared memory uses an atomic xor operation. Each thread in a warp accesses only one bank, avoiding bank conflicts. A p-reduction operation on shared memory is required to combine the 32 partial results.

        
        

-       The maximum number of rows per slice is calculated as \emph{size of shared memory per SM in bits / (32 * blocking factor)} where 32 is the number of threads in a warp. This translates to 192, 96, and 48 rows per slice for blocking factor of 64, 128, and 256 bit, respectively. Hence, one byte per row index is sufficient for this format.       
+       The maximum number of rows per slice is calculated as \emph{size of shared memory per SM in bits} / \emph{(32 * blocking factor)} where 32 is the number of threads in a warp. This translates to 192, 96, and 48 rows per slice for blocking factors of 64, 128, and 256-bits, respectively. Hence, one byte per row index is sufficient for this format.      

        
        

-       \subsubsection*{Large Sliced (LS) COO}
-       In this subformat, the result for each row gets one entry in shared memory, which is shared among all threads in the thread block. Access to shared memory uses an atomic XOR operation. Thus, there will be bank conflicts. However, this drawback can be compensated by a higher texture cache hit rate. There is no p-reduction on shared memory required.
+       \subsubsection*{Large sliced (LS) COO}
+       In this subformat, the result for each row gets one entry in shared memory, which is shared among all threads in the thread block. Access to shared memory uses an atomic xor operation. Thus, there will be bank conflicts. However, this drawback can be compensated for by a higher texture cache hit rate. There is no p-reduction on shared memory required.
+
+       The maximum number of rows per slice is calculated as \emph{size of shared memory per SM in bits} / \emph{blocking factor}. This translates to 6144, 3072, and 1536 rows per slice for blocking factors of 64, 128, and 256-bits, respectively. We need two bytes for the row index.

 
 

-       The maximum number of rows per slice is calculated as \emph{size of shared memory per SM in bits / blocking factor}. This translates to 6144, 3072, and 1536 rows per slice for blocking factor of 64, 128, and 256 bit, respectively. We use two bytes for the row index.
-       \linebreak
-       \linebreak

        \noindent Table \ref{table:scoo} summarizes the differences between the three subformats.
        \begin{itemize}
        \item[a)] Each thread in a block is allocated one memory entry for the result of a row and works on consecutive 12 rows. A p-reduction to Thread\#1 is needed to get the final result of one row.
        \noindent Table \ref{table:scoo} summarizes the differences between the three subformats.
        \begin{itemize}
        \item[a)] Each thread in a block is allocated one memory entry for the result of a row and works on consecutive 12 rows. A p-reduction to Thread\#1 is needed to get the final result of one row.

-       \item[b)] Each warp in a block is allocated one memory entry per row and works on 192 consecutive rows, threads in each warp use atomic XOR to update the value in a given entry. A p-reduction for each row to memory in Warp\#1 is needed to obtain the final result of that row.
-       \item[c)] A thread block works on 6144 consecutive rows and share the memory entry of each row. Any update to the memory has to use atomic XOR.
+       \item[b)] Each warp in a block is allocated one memory entry per row and works on 192 consecutive rows, threads in each warp use atomic xor to update the value in a given entry. A p-reduction for each row to memory in Warp\#1 is needed to obtain the final result of that row.
+       \item[c)] A thread block works on 6144 consecutive rows and shares the memory entry of each row. Any update to the memory has to use atomic xor.

        \end{itemize}
 
        \subsection{Determining the cut-off point of each format}
        \label{ch19:mat-partition-1}
        \end{itemize}
 
        \subsection{Determining the cut-off point of each format}
        \label{ch19:mat-partition-1}

-       To determine which format to use, we compare the performance of two consecutive formats in terms of giga non-zeros (\emph{gnnz}) per second for a given matrix, starting with the dense format and the SS-COO format. The two formats start from the same row (starting from the first row) and work on the minimum number of rows possible. For the dense format, the minimum number of rows is 32. For the SS-COO format (and its variants), the minimum number of rows is \emph{the number of rows in a slice} times \emph{the number of multiprocessors in the GPU}, since one thread block works on one slice and one thread block is assigned to one multiprocessor. 
+       To determine which format to use, we compare the performance of two consecutive formats in terms of giga nonzeros (\emph{gnnz}) per second for a given matrix, starting with the dense format and the SS-COO format. The two formats start from the same row (the first row) and work on the minimum number of rows possible. For the dense format, the minimum number of rows is 32. For the SS-COO format (and its variants), the minimum number of rows is \emph{the number of rows in a slice} times \emph{the number of multiprocessors in the GPU}, since one thread block works on one slice and one thread block is assigned to one multiprocessor. 

        
        The next comparison depends on the result of the current comparison. If the dense format performs better, we decide to use it for rows 1 to 32, and we continue comparing the dense format and the SS-COO format starting from row 33. However, if the SS-COO format performs better, we compare its performance with the next format, MS-COO, starting from the same row, and so on. The idea is to stop considering the denser format once the sparser format outperforms it. Once we get to the comparison between MS-COO and LS-COO, and LS-COO performs better, we don't need to do any further comparisons. LS-COO should be used for the rest of the matrix.
        
        
        The next comparison depends on the result of the current comparison. If the dense format performs better, we decide to use it for rows 1 to 32, and we continue comparing the dense format and the SS-COO format starting from row 33. However, if the SS-COO format performs better, we compare its performance with the next format, MS-COO, starting from the same row, and so on. The idea is to stop considering the denser format once the sparser format outperforms it. Once we get to the comparison between MS-COO and LS-COO, and LS-COO performs better, we don't need to do any further comparisons. LS-COO should be used for the rest of the matrix.
        

-       For Sliced COO format, it is essential to note that when one slice is assigned to each multiprocessor, the load for one multiprocessor may be much higher than that of the other multiprocessors. This is because the matrix rows have been reordered by their weight in a non-increasing order, so the first slice contains more non-zero entries than the rest. Thus, we need to further reorder the rows such that each multiprocessor gets the same level of load.
+       For Sliced COO format, it is essential to note that when one slice is assigned to each multiprocessor, the load for one multiprocessor may be much higher than that of the other multiprocessors. This is because the matrix rows have been reordered by their weight in a nonincreasing order, so the first slice contains more nonzero entries than the rest. Thus, we need to further reorder the rows so that each multiprocessor gets the same level of load.

        
        The cut-off point determination is performed once per matrix and as a preprocessing step. Its runtime is equivalent to less than five iterations of SpMV.
        
        
        
        The cut-off point determination is performed once per matrix and as a preprocessing step. Its runtime is equivalent to less than five iterations of SpMV.
        
        

-\section{SCOO for Single-Precision Floating-Point Matrices}
-The design of SCOO depends on performing atomic operations to combine the partia results. Atomic operations for double-precision are not available on current CUDA-enabled devices. Thus, we only present the SCOO extension for single-precision floating-point matrices.
+\section{SCOO for single-precision floating-point matrices}
+The design of SCOO depends on performing atomic operations to combine the partial results. Atomic operations for double-precision are not available on current CUDA-enabled devices. Thus, we present the SCOO extension for only single-precision floating-point matrices.

 
 

-To extend SCOO for genereal matrices in real numbers, we first adapt to the variety of sparsity in real number matrices. Unlike NFS matrices, the sparsity is not predictable. Hence instead of using three subformats (SS-COO, MS-COO, LS-COO) with a fixed slice size (12, 192, 6144), we make the slice size a variable. Let $h$ denote the slice size, $S$ denote the size in byte of shared memory per thread block ($S=49152$ for Fermi) and $b$ denote the size of each matrix value in byte ($b=4$ for single-precision floating point). $h$ must be a positive integer value of the form of $h = \frac{S}{2^{i} b}$, $\forall{0 \leq i \leq T}$ where $T=\log_2{\textit{(maximum number of thread per block)}}$. 
+To extend SCOO for general matrices in real numbers, we first adapt to the variety of sparsity in real number matrices. Unlike NFS matrices, the sparsity is not predictable. Hence, instead of using three subformats (SS-COO, MS-COO, LS-COO) with a fixed slice size (12, 192, 6144), we make the slice size a variable. Let $h$ denote the slice size, $S$ denote the size in bytes of shared memory per thread block ($S=49152$ for Fermi), and $b$ denote the size of each matrix value in bytes ($b=4$ for single-precision floating point). $h$ must be a positive integer value of the form of $h = \frac{S}{2^{i} b}$, $\forall{0 \leq i \leq T}$ where $T=\log_2{\textit{(maximum number of thread per block)}}$. 

 
 

-For example in Fermi GPUs, single-precision floating point, $h \in \{12288, 6144, 3072, 1536, 768, 384, 192, 96, 48, 24, 12\}$ for $S=49152, b=4, T=10$. The algorithm in Section \ref{ch19:mat-partition-1} can be used to determine the cut-off point for each sub-formats. Figure \ref{fig:final-scoo} illustrates the final matrix in SCOO format assuming the number of multiprocessors in the GPU is $5$ and the row index is ranging from $0..d-1$. $H_i$ and $r_i$ is the slice size and starting row index of the sub-matrix $i$.
+For example in Fermi GPUs, single-precision floating point, $h \in \{12288, 6144, 3072, 1536, 768, 384, 192, 96, 48, 24, 12\}$ for $S=49152, b=4, T=10$. The algorithm described in Section \ref{ch19:mat-partition-1} can be used to determine the cut-off point for each subformat. Figure \ref{fig:final-scoo} illustrates the final matrix in SCOO format assuming the number of multiprocessors in the GPU is $5$ and the row index  ranges from $0..d-1$. $H_i$ and $r_i$ are the slice size and starting row index of the sub-matrix $i$.

 
 \begin{figure}[t]
        \centering
 
 \begin{figure}[t]
        \centering


@@ -279,14 +278,13 @@ For example in Fermi GPUs, single-precision floating point, $h \in \{12288, 6144
        \label{fig:final-scoo}
 \end{figure}
 
        \label{fig:final-scoo}
 \end{figure}
 

-More details of the SCOO format for floating-point matrices can be found in \cite{ch19:spmv-iccs}
+More details of the SCOO format for floating-point matrices can be found in \cite{ch19:spmv-iccs}.

 
 

-\section{Performance Evaluation}
+\section{Performance evaluation}

 \label{ch19:result}
 
 \label{ch19:result}
 

-\subsection{Experimental Setup}
+\subsection{Experimental setup}

 \begin{table}
 \begin{table}

-\caption{Overview of hardware used in the experiments}

 \begin{center}
 \begin{tabular}{|l|l|l|l|}
 \hline
 \begin{center}
 \begin{tabular}{|l|l|l|l|}
 \hline


@@ -304,21 +302,20 @@ Max Memory bandwidth (GB/s)& 144 & 192 & 21 \\
 \hline
 \end{tabular}
 \end{center}
 \hline
 \end{tabular}
 \end{center}

+\caption{Overview of hardware used in the experiments}

 \label{table:hardware}
 \end{table} 
 
 \label{table:hardware}
 \end{table} 
 

-Table \ref{table:hardware} specifies the GPUs and the CPU workstation used for performance evaluation. The performance is measured in terms of $ggnz$ and $Gflop/s$ for binary and floating-point matrices. Measured GPU performance does neither include PCIe data transfers nor matrix preprocessing. These are realistic assumptions since SpMV applications usually consist of a large number of iterations where the sparse matrix is iteratively multiplied by the input/output vectors.
+Table \ref{table:hardware} specifies the GPUs and the CPU workstation used for performance evaluation. The performance is measured in terms of $gnnz$ and $Gflop/s$ for binary and floating-point matrices. Measured GPU performance includes neither PCIe data transfers nor matrix preprocessing. These are realistic assumptions since SpMV applications usually consist of a large number of iterations where the sparse matrix is iteratively multiplied by the input/output vectors.

 
 Our source code is compiled and executed using CUDA toolkit and driver version $4.2$ under Linux Ubuntu 12.04.
 
 \subsection{Experimental results of SpMV on NFS matrix}
 We have evaluated our implementation on an NVIDIA Tesla C2075 (ECC disabled) with 6 GB RAM. We have compared the GPU performance with the open-source CADO-NFS \cite{ch19:cadonfs} program running on Intel Core i7-920 CPU with 12 GB DDR3-1066 memory. The RSA-170 matrix (see Table \ref{table:rsa_matrix}) is used for performance evaluation. 
 
 
 Our source code is compiled and executed using CUDA toolkit and driver version $4.2$ under Linux Ubuntu 12.04.
 
 \subsection{Experimental results of SpMV on NFS matrix}
 We have evaluated our implementation on an NVIDIA Tesla C2075 (ECC disabled) with 6 GB RAM. We have compared the GPU performance with the open-source CADO-NFS \cite{ch19:cadonfs} program running on Intel Core i7-920 CPU with 12 GB DDR3-1066 memory. The RSA-170 matrix (see Table \ref{table:rsa_matrix}) is used for performance evaluation. 
 

-The speedups compared to the multi-threaded CADO-NFS bucket implementation on an Intel Core i7-920 are given in brackets. The GPU memory required to store the sparse matrix and the corresponding bytes per non-zero (nnz) are also reported. CADO-NFS contains several CPU optimized SpMV implementations using multi-threading and SSE instructions. In this experiment, we compare the performance of our implementation to the CPU cache-optimized bucket format of CADO-NFS using 8 threads.
+The speedups compared to the multithreaded CADO-NFS bucket implementation on an Intel Core i7-920 are given in parentheses in Table~\ref{table:hardware}. The GPU memory required to store the sparse matrix and the corresponding bytes per nonzero (nnz) are also reported. CADO-NFS contains several CPU optimized SpMV implementations using multithreading and SSE instructions. In this experiment, we compare the performance of our implementation to the CPU cache-optimized bucket format of CADO-NFS using 8 threads.

 
 \begin{table}[t]
 
 \begin{table}[t]

-       \caption{Performance of SpMV on RSA-170 matrix}
-       \label{table:rsa170}

        \centering
        \begin{tabular}
        { | c || c | c | c | c |}       
        \centering
        \begin{tabular}
        { | c || c | c | c | c |}       


@@ -332,17 +329,17 @@ The speedups compared to the multi-threaded CADO-NFS bucket implementation on an
 
        \hline
        \end{tabular}
 
        \hline
        \end{tabular}

+       \caption{Performance of SpMV on RSA-170 matrix}
+       \label{table:rsa170}

 \end{table}
                        
        Table \ref{table:rsa170} shows the result. The dual-GPU implementation achieves speedups between 1.93 and 1.96 compared to the single-GPU performance.
        
 \end{table}
                        
        Table \ref{table:rsa170} shows the result. The dual-GPU implementation achieves speedups between 1.93 and 1.96 compared to the single-GPU performance.
        

-       Table \ref{table:individual} show the individual performance of each sub-format (in \emph{gnnz/s}), the percentage of non-zeros are included in bracket. The result show that the performance decreases when the matrix gets sparser. The MS-COO and LS-COO performance degrade when the blocking factor is increased from 64 to 128 and 256 bit. This is caused by the increased number of bank conflicts and serialization of atomic XOR operations on larger blocking factors. Thus, the SS-COO format gets a higher percentage of non-zeros with 128 and 256 bit blocking factor. 
+       Table \ref{table:individual} shows the individual performance of each subformat (in \emph{gnnz/s}), the percentage of nonzeros are included in parentheses. The results show that the performance decreases when the matrix gets sparser. The MS-COO and LS-COO performance degrade when the blocking factor is increased from 64 to 128 and 256-bit. This is caused by the increased number of bank conflicts and serialization of atomic xor operations on larger blocking factors. Thus, the SS-COO format gets a higher percentage of nonzeros with 128 and 256-bit blocking factor. 

        
        

-       More detail results of our full Block Wiedemann CUDA implementation as well as a multi-GPU implementation can be found in \cite{ch19:spmv-ccpe}
+       More detail results of our full block Wiedemann CUDA implementation as well as a multi-GPU implementation can be found in \cite{ch19:spmv-ccpe}

 
 \begin{table}[h]
 
 \begin{table}[h]

-       \caption{Performance for each of the four sub-format partitions of the RSA-170 matrix on a C2075.}
-       \label{table:individual}

        \centering
        \begin{tabular}
        { | c || c | c | c | c | c |}   
        \centering
        \begin{tabular}
        { | c || c | c | c | c | c |}   


@@ -351,17 +348,27 @@ The speedups compared to the multi-threaded CADO-NFS bucket implementation on an
        \hline            
          64     & 13.66 (24\%) & 9.66 (11\%) & 8.13 (13\%) & 2.77 (52\%)  \\
          128    & 9.66 (15\%) & 7.53 (24\%) & 3.23 (6\%) & 1.86 (55\%) \\
        \hline            
          64     & 13.66 (24\%) & 9.66 (11\%) & 8.13 (13\%) & 2.77 (52\%)  \\
          128    & 9.66 (15\%) & 7.53 (24\%) & 3.23 (6\%) & 1.86 (55\%) \\

-         256    & 7.00 (15\%) & 4.21 (21\%) & 2.34 (6\%) & 1.33 (58 \%) \\               
+         256    & 7.00 (15\%) & 4.21 (21\%) & 2.34 (6\%) & 1.33 (58\%) \\                

        \hline
        \end{tabular}
        \hline
        \end{tabular}

+       \caption{Performance for each of the four subformat partitions of the RSA-170 matrix on a C2075.}
+       \label{table:individual}

 \end{table}
 
 \subsection{Experimental results of SpMV on floating-point matrices}
 \label{ch19:result-2}
 \npthousandsep{,}
 
 \end{table}
 
 \subsection{Experimental results of SpMV on floating-point matrices}
 \label{ch19:result-2}
 \npthousandsep{,}
 

+
+
+
+
+\subsubsection{Performance comparison to existing GPU formats}
+We compare the performance of our SCOO format to available SpMV implementations on both GPU and CPU. The set of selected test matrices are collected from the University of Florida Sparse Matrix Collection \cite{ch19:matrix-collection}. We have chosen the biggest matrices from different areas that with their corresponding input and output vectors can still fit into the 6GB global memory of a C2075 GPU. Table \ref{table:matrices} gives an overview of those matrices. 
+
+
+

 \begin{table}
 \begin{table}

-\caption{Overview of sparse matrices used for performance evaluation}

 \begin{center}
 \begin{tabular}{|l|n{8}{0}|n{9}{0}|n{3}{2}|l|}
 \hline
 \begin{center}
 \begin{tabular}{|l|n{8}{0}|n{9}{0}|n{3}{2}|l|}
 \hline


@@ -397,11 +404,32 @@ uk-2002 &18520486&18520486&16.10&directed graph\\
 \hline
 \end{tabular}
 \end{center}
 \hline
 \end{tabular}
 \end{center}

+\caption{Overview of sparse matrices used for performance evaluation.}

 \label{table:matrices}
 \end{table}
 
 \label{table:matrices}
 \end{table}
 

-\begin{figure}
-\label{fig:mat-visual}
+We compare the SCOO format to the CSR, COO, and HYB format of Cusp 0.3.0. Other Cusp formats are not able to run on the large tested matrices that we selected. The results are shown in Figure \ref{fig:scoo-vs-gpu}. The performances are in terms of Gflop/s which is based on the assumption of two flops per nonzero entry of the matrix \cite{ch19:nvidia-spmv, ch19:bellpack}.
+
+
+
+\begin{figure}[htbp]
+\centering
+\includegraphics[height=5cm]{Chapters/chapter19/fig/scoo-a.pdf}
+\caption{Performance comparison of SCOO and other GPU formats for each test matrix on a Fermi Tesla C2075 (ECC disabled).}
+\label{fig:scoo-vs-gpu}
+\end{figure}
+
+
+
+The SCOO format achieves a stable performance for different matrices in single-precision mode. In most cases a performance of over 10 Gflop/s can be sustained. For some highly unstructured matrices such as \emph{GL7d19}, \emph{wikipedia-20070206}, \emph{rgg\_n\_2\_24\_s0} and \emph{kron\_g500-logn21} SCOO achieves high speedups ranging from 3 to 6 compared to the best performaning Cusp format. 
+
+For most matrices, HYB produces the best performance among the tested Cusp formats. HYB is able to outperform SCOO only for two matrices: \emph{nlpkkt120} and \emph{nlpkkt160}. Both matrices have a similar structure i.e. they consist of consecutive rows that have a very similar number of nonzero coefficients which is suitable to be stored in the ELL section of the HYB format. Moreover the nonzeros are close to each other facilitating coaleasing and cache-friendly access patterns by nature. SCOO is able to outperform COO and CSR for all tested matrices.
+
+In matrix $Relat9$ we observe some patterns but the matrix is still generally unstructured, thus SCOO is able to achieve about 2 times speed up compared to HYB which is the best among tested Cusp formats in this case. The average speedup of SCOO for the tested matrices is 3.0 compared to CSR, 5.02 compared to COO, 2.15 compared to HYB.
+
+We show the visualization of sparse matrices \emph{nlpkkt120}, \emph{Relat9}, \emph{GL7d19} in Figure \ref{fig:mat-str}, \ref{fig:mat-mid}, \ref{fig:mat-unstr} using MatView \cite{ch19:matview}. The white color represents zero entries, gray color represents nonzero entries.
+
+\begin{figure}[htbp]

 \centering
        \subfigure[nlpkkt120] {
                \includegraphics[width=100pt]{Chapters/chapter19/fig/matrix-str.pdf}
 \centering
        \subfigure[nlpkkt120] {
                \includegraphics[width=100pt]{Chapters/chapter19/fig/matrix-str.pdf}


@@ -417,32 +445,14 @@ uk-2002 &18520486&18520486&16.10&directed graph\\
        }
   
   \caption{Visualization of \emph{nlpkkt120}, \emph{Relat9}, and \emph{GL7d19} matrix.}
        }
   
   \caption{Visualization of \emph{nlpkkt120}, \emph{Relat9}, and \emph{GL7d19} matrix.}

+%  \label{fig:mat-visual}

 \end{figure}
 
 \end{figure}
 

-\begin{figure}
-\centering
-\includegraphics[height=5cm]{Chapters/chapter19/fig/scoo-a.pdf}
-\caption{Performance comparison of SCOO and other GPU formats for each test matrix on a Fermi Tesla C2075 (ECC disabled).}
-\label{fig:scoo-vs-gpu}
-\end{figure}
-
-\subsubsection{Performance comparison to existing GPU formats}
-We compare the performance of our SCOO format to available SpMV implementations on both GPU and CPU. The set of selected test matrices are collected from the University of Florida Sparse Matrix Collection \cite{ch19:matrix-collection}. We have chosen the biggest matrices from different areas that with their corresponding input and output vector can still fit into the 6GB global memory of a C2075 GPU. Table \ref{table:matrices} gives an overview of those matrices. 
-
-We compare the SCOO format to the CSR, COO and HYB format of Cusp 0.3.0. Other Cusp formats are not able to run on the large tested matrices that we selected. The results are shown in Figure \ref{fig:scoo-vs-gpu}. The performance are in terms of Gflop/s which is based on the assumption of two flops per non-zero entry of the matrix \cite{ch19:nvidia-spmv, ch19:bellpack}.
-
-The SCOO format achieves a stable performance for different matrices in single-precision mode. In most cases a performance of over 10 Gflop/s can be sustained. For some highly unstructured matrices such as \emph{GL7d19}, \emph{wikipedia-20070206}, \emph{rgg\_n\_2\_24\_s0} and \emph{kron\_g500-logn21} SCOO achieves high speedups ranging from 3 to 6 compared to the best performaning Cusp format. 
-
-For most matrices, HYB produces the best performance among the tested Cusp formats. HYB is able to outperform SCOO only for two matrices: \emph{nlpkkt120} and \emph{nlpkkt160}. Both matrices have a similar structure i.e. they consist of consecutive rows that have a very similar number of non-zero coefficients which is suitable to be stored in the ELL section of the HYB format. Moreover the non-zeros are close to each other facilitating coaleasing and cache-friendly access patterns by nature. SCOO is able to outperform COO and CSR for all tested matrices.
-
-In matrix $Relat9$ we observe some patterns but the matrix is still generally unstructured, thus SCOO is able to achieve about 2 times speed up compared to HYB which is the best among tested Cusp formats in this case. The average speedup of SCOO for the tested matrices is 3.0 compared to CSR, 5.02 compared to COO, 2.15 compared to HYB.
-
-We show the visualization of sparse matrices \emph{nlpkkt120}, \emph{Relat9}, \emph{GL7d19} in Figure \ref{fig:mat-str}, \ref{fig:mat-mid}, {fig:mat-unstr} using MatView \cite{ch19:matview}. The white color represents zero entries, gray color represents non-zero entries.

 
 \subsubsection{Performance comparison to a CPU implementation}
 \begin{figure}
 \centering
 
 \subsubsection{Performance comparison to a CPU implementation}
 \begin{figure}
 \centering

-\includegraphics[width=200pt]{Chapters/chapter19/fig/gpu-cpu.pdf}
+\includegraphics[width=340pt]{Chapters/chapter19/fig/gpu-cpu.pdf}

 
 \caption{Performance of the SCOO on a GTX-580 and a CPU implementation using MKL performed on a Core-i7 2700K using 8 threads.}
 \label{fig:scoo-vs-cpu}
 
 \caption{Performance of the SCOO on a GTX-580 and a CPU implementation using MKL performed on a Core-i7 2700K using 8 threads.}
 \label{fig:scoo-vs-cpu}

diff --git a/BookGPU/Chapters/chapter3/ch3.aux b/BookGPU/Chapters/chapter3/ch3.aux
index 0fb9a0d1efe6382a59824bb7541d222844791f8d..8740459458f83a1369dab80e44d6feba75519050 100644 (file)
--- a/BookGPU/Chapters/chapter3/ch3.aux
+++ b/BookGPU/Chapters/chapter3/ch3.aux
@@ -43,14 +43,6 @@
 \@writefile{lot}{\contentsline {table}{\numberline {4.1}{\ignorespaces Performance results of \texttt  {kernel medianR}. \relax }}{35}}
 \newlabel{tab:medianHisto1}{{4.1}{35}}
 \@writefile{toc}{\contentsline {section}{\numberline {4.3}NVIDIA GPU tuning recipes}{35}}
 \@writefile{lot}{\contentsline {table}{\numberline {4.1}{\ignorespaces Performance results of \texttt  {kernel medianR}. \relax }}{35}}
 \newlabel{tab:medianHisto1}{{4.1}{35}}
 \@writefile{toc}{\contentsline {section}{\numberline {4.3}NVIDIA GPU tuning recipes}{35}}

-\newlabel{img:sap_example_ref}{{4.3(a)}{36}}
-\newlabel{sub@img:sap_example_ref}{{(a)}{36}}
-\newlabel{img:sap_example_med3}{{4.3(b)}{36}}
-\newlabel{sub@img:sap_example_med3}{{(b)}{36}}
-\newlabel{img:sap_example_med5}{{4.3(c)}{36}}
-\newlabel{sub@img:sap_example_med5}{{(c)}{36}}
-\newlabel{img:sap_example_med3_it2}{{4.3(d)}{36}}
-\newlabel{sub@img:sap_example_med3_it2}{{(d)}{36}}

 \@writefile{lof}{\contentsline {figure}{\numberline {4.3}{\ignorespaces Example of median filtering, applied to salt and pepper noise reduction.\relax }}{36}}
 \@writefile{lof}{\contentsline {subfigure}{\numberline{(a)}{\ignorespaces {Airplane image, corrupted by salt and pepper noise of density 0.25}}}{36}}
 \@writefile{lof}{\contentsline {subfigure}{\numberline{(b)}{\ignorespaces {Image denoised by a $3\times 3$ median filter}}}{36}}
 \@writefile{lof}{\contentsline {figure}{\numberline {4.3}{\ignorespaces Example of median filtering, applied to salt and pepper noise reduction.\relax }}{36}}
 \@writefile{lof}{\contentsline {subfigure}{\numberline{(a)}{\ignorespaces {Airplane image, corrupted by salt and pepper noise of density 0.25}}}{36}}
 \@writefile{lof}{\contentsline {subfigure}{\numberline{(b)}{\ignorespaces {Image denoised by a $3\times 3$ median filter}}}{36}}


@@ -92,14 +84,6 @@
 \@writefile{toc}{\contentsline {subsection}{\numberline {4.5.2}Fast approximated $n\times n$ median filter }{47}}
 \@writefile{lot}{\contentsline {table}{\numberline {4.3}{\ignorespaces Measured performance of one generic pseudo-separable median kernel applied to 4096$\times $4096 pixel image with various window sizes.\relax }}{48}}
 \newlabel{tab:medianSeparable}{{4.3}{48}}
 \@writefile{toc}{\contentsline {subsection}{\numberline {4.5.2}Fast approximated $n\times n$ median filter }{47}}
 \@writefile{lot}{\contentsline {table}{\numberline {4.3}{\ignorespaces Measured performance of one generic pseudo-separable median kernel applied to 4096$\times $4096 pixel image with various window sizes.\relax }}{48}}
 \newlabel{tab:medianSeparable}{{4.3}{48}}

-\newlabel{img:sap_example_ref}{{4.11(a)}{49}}
-\newlabel{sub@img:sap_example_ref}{{(a)}{49}}
-\newlabel{img:sap_example_sep_med3}{{4.11(b)}{49}}
-\newlabel{sub@img:sap_example_sep_med3}{{(b)}{49}}
-\newlabel{img:sap_example_sep_med5}{{4.11(c)}{49}}
-\newlabel{sub@img:sap_example_sep_med5}{{(c)}{49}}
-\newlabel{img:sap_example_sep_med3_it2}{{4.11(d)}{49}}
-\newlabel{sub@img:sap_example_sep_med3_it2}{{(d)}{49}}

 \@writefile{lof}{\contentsline {figure}{\numberline {4.11}{\ignorespaces Example of separable median filtering (smoother), applied to salt and pepper noise reduction.\relax }}{49}}
 \@writefile{lof}{\contentsline {subfigure}{\numberline{(a)}{\ignorespaces {Airplane image, corrupted with by salt and pepper noise of density 0.25}}}{49}}
 \@writefile{lof}{\contentsline {subfigure}{\numberline{(b)}{\ignorespaces {Image denoised by a $3\times 3$ separable smoother}}}{49}}
 \@writefile{lof}{\contentsline {figure}{\numberline {4.11}{\ignorespaces Example of separable median filtering (smoother), applied to salt and pepper noise reduction.\relax }}{49}}
 \@writefile{lof}{\contentsline {subfigure}{\numberline{(a)}{\ignorespaces {Airplane image, corrupted with by salt and pepper noise of density 0.25}}}{49}}
 \@writefile{lof}{\contentsline {subfigure}{\numberline{(b)}{\ignorespaces {Image denoised by a $3\times 3$ separable smoother}}}{49}}

diff --git a/BookGPU/Chapters/chapter3/ch3.tex b/BookGPU/Chapters/chapter3/ch3.tex
index 64709ea4e51f510d3925a37e43f552f5852e8068..7078aaffe378b04beaeef5bbda0dfe0aa55407f5 100755 (executable)
--- a/BookGPU/Chapters/chapter3/ch3.tex
+++ b/BookGPU/Chapters/chapter3/ch3.tex
@@ -181,10 +181,10 @@ On the GPU's side, we note high dependence on window size due to the redundancy
 
 \begin{figure}[t]
 \centering
 
 \begin{figure}[t]
 \centering

-   \subfigure[Airplane image, corrupted by salt and pepper noise of density 0.25]{\label{img:sap_example_ref} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25.png}}\qquad
-   \subfigure[Image denoised by a $3\times 3$ median filter]{\label{img:sap_example_med3} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_med3.png}}\\
-   \subfigure[Image denoised by a $5\times 5$ median filter]{\label{img:sap_example_med5} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_med5.png}}\qquad
-   \subfigure[Image denoised by 2 iterations of a $3\times 3$ median filter]{\label{img:sap_example_med3_it2} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_med3_it2.png}}\\
+   \subfigure[Airplane image, corrupted by salt and pepper noise of density 0.25]{ \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25.png}}\qquad
+   \subfigure[Image denoised by a $3\times 3$ median filter]{ \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_med3.png}}\\
+   \subfigure[Image denoised by a $5\times 5$ median filter]{ \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_med5.png}}\qquad
+   \subfigure[Image denoised by 2 iterations of a $3\times 3$ median filter]{ \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_med3_it2.png}}\\

    \caption{Example of median filtering, applied to salt and pepper noise reduction.}
    \label{fig:sap_examples}
 \end{figure}
    \caption{Example of median filtering, applied to salt and pepper noise reduction.}
    \label{fig:sap_examples}
 \end{figure}


@@ -351,10 +351,10 @@ Torben Morgensen sorting takes place between lines 37 and 66 and eventually, the
 It has to be noticed that this smoother, unlike the technique we proposed for fixed-size median filters, cannot be considered as a state-of-the-art technique as, for example, the one presented in \cite{4287006}. However, it may be considered as a good, easy to use and efficient alternative as confirmed by the results presented in Table \ref{tab:medianSeparable}. Pixel throughput values achieved by our kernel, though not constant with window size, remain very competitive if window size is kept under $120\times 120$ pixels, especially when outputting 2 pixels per thread (in \cite{4287006}, pixel throughput is around 7MP/s).
 Figure \ref{fig:sap_examples2} shows an example of a $512\times 512$ pixel image, corrupted by a  \textit{salt and pepper} noise, and the denoised versions, outputted respectively by a $3\times 3$, a $5\times 5$, and a $55\times 55 $ separable smoother.
 \begin{figure}
 It has to be noticed that this smoother, unlike the technique we proposed for fixed-size median filters, cannot be considered as a state-of-the-art technique as, for example, the one presented in \cite{4287006}. However, it may be considered as a good, easy to use and efficient alternative as confirmed by the results presented in Table \ref{tab:medianSeparable}. Pixel throughput values achieved by our kernel, though not constant with window size, remain very competitive if window size is kept under $120\times 120$ pixels, especially when outputting 2 pixels per thread (in \cite{4287006}, pixel throughput is around 7MP/s).
 Figure \ref{fig:sap_examples2} shows an example of a $512\times 512$ pixel image, corrupted by a  \textit{salt and pepper} noise, and the denoised versions, outputted respectively by a $3\times 3$, a $5\times 5$, and a $55\times 55 $ separable smoother.
 \begin{figure}

-   \subfigure[Airplane image, corrupted with by salt and pepper noise of density 0.25]{\label{img:sap_example_ref} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25.png}}\qquad
-   \subfigure[Image denoised by a $3\times 3$ separable smoother]{\label{img:sap_example_sep_med3} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_sep_med3.png}}\\
-   \subfigure[Image denoised by a $5\times 5$ separable smoother]{\label{img:sap_example_sep_med5} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_sep_med5.png}}\qquad
-   \subfigure[Image background estimation by a $55\times 55$ separable smoother]{\label{img:sap_example_sep_med3_it2} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_sep_med111.png}}\\
+   \subfigure[Airplane image, corrupted with by salt and pepper noise of density 0.25]{ \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25.png}}\qquad
+   \subfigure[Image denoised by a $3\times 3$ separable smoother]{ \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_sep_med3.png}}\\
+   \subfigure[Image denoised by a $5\times 5$ separable smoother]{ \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_sep_med5.png}}\qquad
+   \subfigure[Image background estimation by a $55\times 55$ separable smoother]{ \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_sep_med111.png}}\\

    \caption{Example of separable median filtering (smoother), applied to salt and pepper noise reduction.}
    \label{fig:sap_examples2}
 \end{figure}
    \caption{Example of separable median filtering (smoother), applied to salt and pepper noise reduction.}
    \label{fig:sap_examples2}
 \end{figure}

diff --git a/BookGPU/Chapters/chapter5/ch5.tex b/BookGPU/Chapters/chapter5/ch5.tex
index e37a709483e7a2e04589219214fcfb4219e10e70..8ee777552706c20728390833626329328a6e87ac 100644 (file)
--- a/BookGPU/Chapters/chapter5/ch5.tex
+++ b/BookGPU/Chapters/chapter5/ch5.tex
@@ -222,8 +222,8 @@ All items are either directly available in the library or can be designed from c
 
 With the right hand side operator in place, we are ready to implement the solver. For this simple PDE problem we compute all necessary initial data in the body of the main function and use the forward Euler time integrator to compute the solution till $t=t_{end}$. For more advanced PDE solvers, a built-in \texttt{ode\_solver} class is defined, that helps taking care of initialization and storage of multiple state variables. Declaring type definitions for all components at the beginning of the main file gives a good overview of the solver composition. In this way, it will be easy to control or change solver components at later times. Listing~\ref{ch5:lst:heattypedefs} lists the type definitions\index{type definitions} that are used to assemble the heat conduction solver. %\stefan{elaborate on the choices}.
 
 
 With the right hand side operator in place, we are ready to implement the solver. For this simple PDE problem we compute all necessary initial data in the body of the main function and use the forward Euler time integrator to compute the solution till $t=t_{end}$. For more advanced PDE solvers, a built-in \texttt{ode\_solver} class is defined, that helps taking care of initialization and storage of multiple state variables. Declaring type definitions for all components at the beginning of the main file gives a good overview of the solver composition. In this way, it will be easy to control or change solver components at later times. Listing~\ref{ch5:lst:heattypedefs} lists the type definitions\index{type definitions} that are used to assemble the heat conduction solver. %\stefan{elaborate on the choices}.
 

-\lstset{label=ch5:lst:heattypedefs,caption={Type definitions for all the heat conduction solver components used throughout the remaining code.}}
-\begin{lstlisting}
+\lstset{caption={Type definitions for all the heat conduction solver components used throughout the remaining code.}}
+\begin{lstlisting}[label=ch5:lst:heattypedefs]

 typedef double                              value_type;
 typedef laplacian<value_type>               rhs_type;
 typedef gpulab::grid<value_type>            vector_type;
 typedef double                              value_type;
 typedef laplacian<value_type>               rhs_type;
 typedef gpulab::grid<value_type>            vector_type;


@@ -232,8 +232,8 @@ typedef gpulab::integration::forward_euler  time_integrator_type;
 \end{lstlisting}
 The grid is by default treated as a device object and memory is allocated on the GPU upon initialization of the grid. Setting up the grid can be done via the property type class. The property class holds information about the discrete and physical dimensions, along with fictitious ghost (halo) layers and periodicity conditions. For the heat conduction problem we use a non-periodic domain within the unit square and with no ghost layers. Listing \ref{ch5:lst:gridsetup} illustrates the grid assembly.
 
 \end{lstlisting}
 The grid is by default treated as a device object and memory is allocated on the GPU upon initialization of the grid. Setting up the grid can be done via the property type class. The property class holds information about the discrete and physical dimensions, along with fictitious ghost (halo) layers and periodicity conditions. For the heat conduction problem we use a non-periodic domain within the unit square and with no ghost layers. Listing \ref{ch5:lst:gridsetup} illustrates the grid assembly.
 

-\lstset{label=ch5:lst:gridsetup,caption={Create 2D grid of discrete dimension \texttt{N} times \texttt{N} and physical dimension $0$ to $1$.}}
-\begin{lstlisting}
+\lstset{caption={Create 2D grid of discrete dimension \texttt{N} times \texttt{N} and physical dimension $0$ to $1$.}}
+\begin{lstlisting}[label=ch5:lst:gridsetup]

 // Setup discrete and physical dimensions
 gpulab::grid_dim<int>          dim(N,N,1);
 gpulab::grid_dim<value_type>   p0(0,0);
 // Setup discrete and physical dimensions
 gpulab::grid_dim<int>          dim(N,N,1);
 gpulab::grid_dim<value_type>   p0(0,0);


@@ -245,8 +245,8 @@ vector_type                    u(props);
 \end{lstlisting}
 Hereafter the vector \texttt{u} can be initialized accordingly to \eqref{ch5:eq:heatinit}. Finally we need to instantiate the right hand side Laplacian operator from Listing \ref{ch5:lst:laplaceimpl} and the forward Euler time integrator in order to integrate from $t_0$ till $t_{end}$.
 
 \end{lstlisting}
 Hereafter the vector \texttt{u} can be initialized accordingly to \eqref{ch5:eq:heatinit}. Finally we need to instantiate the right hand side Laplacian operator from Listing \ref{ch5:lst:laplaceimpl} and the forward Euler time integrator in order to integrate from $t_0$ till $t_{end}$.
 

-\lstset{label=ch5:lst:timeintegrator,caption=Create time integrator and the right hand side Laplacian operator.}
-\begin{lstlisting}
+\lstset{caption=Create time integrator and the right hand side Laplacian operator.}
+\begin{lstlisting}[label=ch5:lst:timeintegrator]

 rhs_type rhs(alpha);           // Create right hand side operator
 time_integrator_type solver;   // Create time integrator
 solver(&rhs,u,0.0f,tend,dt);   // Integrate from to tend using dt
 rhs_type rhs(alpha);           // Create right hand side operator
 time_integrator_type solver;   // Create time integrator
 solver(&rhs,u,0.0f,tend,dt);   // Integrate from to tend using dt


@@ -304,8 +304,8 @@ where $k$ is the iteration number and $\mymat{M}$ is the preconditioner\index{pr
 
 In the following section we demonstrate how to assemble a solver for the discrete Poisson problem, using one of the tree iterative methods to efficiently solve \eqref{ch5:eq:poissonsystem}.
 
 
 In the following section we demonstrate how to assemble a solver for the discrete Poisson problem, using one of the tree iterative methods to efficiently solve \eqref{ch5:eq:poissonsystem}.
 

-\lstset{label=ch5:lst:dc,caption={Main loop for the iterative defect correction solver. The solver is instantiated with template argument types for the matrix and vector classes, allowing underlying implementations to be based on GPU kernels.}}
-\begin{lstlisting}
+\lstset{caption={Main loop for the iterative defect correction solver. The solver is instantiated with template argument types for the matrix and vector classes, allowing underlying implementations to be based on GPU kernels.}}
+\begin{lstlisting}[label=ch5:lst:dc]

 while(r.nrm2() > tol)
 {
     // Calculate residual
 while(r.nrm2() > tol)
 {
     // Calculate residual


@@ -330,8 +330,8 @@ Assembling the Poisson solver follows almost the same procedure as the heat cond
 
 At the beginning of the solver implementation we list the type definitions for the Poisson solver that will be used throughout the implementation. Here we use a geometric multigrid\index{multigrid} method as a preconditioner for the defect correction method. Therefore the multigrid solver is assembled first, so that it can be used in the assembling of the defect correction method. Listing \ref{ch5:lst:poissontypedefs} defines the types for the vector, the matrix, the multigrid preconditioner and the defect correction solver. The geometric multigrid method needs two additional template arguments, that are specific for multigrid, namely a smoother and a grid restriction/interpolation operator. These arguments are free to be implemented and supplied by the developer if special care is required for their specific problems, e.g. for a custom grid structure. For the Poisson problem on a regular grid, the library contains built-in restriction and interpolation operators, and a red-black Gauss-Seidel smoother. We refer to \cite{ch5:Trottenberg2001} for extensive details on multigrid methods. The monitor and config types that appear in Listing \ref{ch5:lst:poissontypedefs} are used for convergence monitoring within the iterative solver and to control run time parameters, such as tolerances, iteration limits, etc.
 
 
 At the beginning of the solver implementation we list the type definitions for the Poisson solver that will be used throughout the implementation. Here we use a geometric multigrid\index{multigrid} method as a preconditioner for the defect correction method. Therefore the multigrid solver is assembled first, so that it can be used in the assembling of the defect correction method. Listing \ref{ch5:lst:poissontypedefs} defines the types for the vector, the matrix, the multigrid preconditioner and the defect correction solver. The geometric multigrid method needs two additional template arguments, that are specific for multigrid, namely a smoother and a grid restriction/interpolation operator. These arguments are free to be implemented and supplied by the developer if special care is required for their specific problems, e.g. for a custom grid structure. For the Poisson problem on a regular grid, the library contains built-in restriction and interpolation operators, and a red-black Gauss-Seidel smoother. We refer to \cite{ch5:Trottenberg2001} for extensive details on multigrid methods. The monitor and config types that appear in Listing \ref{ch5:lst:poissontypedefs} are used for convergence monitoring within the iterative solver and to control run time parameters, such as tolerances, iteration limits, etc.
 

-\lstset{label=ch5:lst:poissontypedefs,caption={Type definitions for the Laplacian matrix component and the multigrid preconditioned iterative defect correction\index{defect correction} solver.}}
-\begin{lstlisting}
+\lstset{caption={Type definitions for the Laplacian matrix component and the multigrid preconditioned iterative defect correction\index{defect correction} solver.}}
+\begin{lstlisting}[label=ch5:lst:poissontypedefs]

 typedef double                                   value_type;
 typedef gpulab::grid<value_type>                 vector_type;
 typedef laplacian<vector_type>                   matrix_type;
 typedef double                                   value_type;
 typedef gpulab::grid<value_type>                 vector_type;
 typedef laplacian<vector_type>                   matrix_type;


@@ -357,8 +357,8 @@ typedef gpulab::solvers::defect_correction<dc_types> dc_solver_type;
 
 With the type definitions set up, the implementation for the Poisson solver follows in Listing \ref{ch5:lst:poissonsolver}. Some of the initializations are left out, as they follow the same procedure as for the Heat conduction example. The defect correction and geometric multigrid solvers are initiated and then the multigrid is set as a preconditioner to the defect correction method. Finally the system is solved via a call to \texttt{solve()}.
 
 
 With the type definitions set up, the implementation for the Poisson solver follows in Listing \ref{ch5:lst:poissonsolver}. Some of the initializations are left out, as they follow the same procedure as for the Heat conduction example. The defect correction and geometric multigrid solvers are initiated and then the multigrid is set as a preconditioner to the defect correction method. Finally the system is solved via a call to \texttt{solve()}.
 

-\lstset{label=ch5:lst:poissonsolver,caption={Initializing the preconditioned defect correction solver to approximate the solution to $\mymat{A}\myvec{u}=\myvec{f}$.}}
-\begin{lstlisting}
+\lstset{caption={Initializing the preconditioned defect correction solver to approximate the solution to $\mymat{A}\myvec{u}=\myvec{f}$.}}
+\begin{lstlisting}[label=ch5:lst:poissonsolver]

 matrix_type A(alpha);               // High-order matrix
 matrix_type M(1);                   // Low-order matrix
 
 matrix_type A(alpha);               // High-order matrix
 matrix_type M(1);                   // Low-order matrix
 


@@ -508,8 +508,8 @@ Using the defined $\mathcal{F}_{\Delta T}$ and $\mathcal{G}_{\Delta T}$ operator
 U_{n}^{k+1}=\mathcal{\mathcal{G}}_{\Delta T}\left(U_{n-1}^{k+1}\right)+\mathcal{\mathcal{F}}_{\Delta T}\left(U_{n-1}^{k}\right)-\mathcal{\mathcal{G}}_{\Delta T}\left(U_{n-1}^{k}\right),\quad U_{0}^{k}=u^{0},
 \end{equation}
 with the initial prediction $U^{0}_{n} = \mathcal{G}_{\Delta T}^{n} u^{0}$ for $n=1\ldots N$ and $k=1\ldots K$. $N$ being the number of time subdomains, while $K\geq1$ is the number of predictor-corrector iterations applied. The parareal algorithm is implemented in the library as a separate time integration component, using a fully distributed work scheduling model, as proposed by Aubanel (2010)~\cite{ch5:EA10}. The model is schematically presented in Figure \ref{ch5:fig:FullyDistributedCores}. The parareal component hides all communication and work distribution from the application developer.  It is defined such that a user only has to decide what coarse and fine propagators to use. Setting up the type definitions for parareal time integration using forward Euler for coarse propagation and fourth order Runge-Kutta for fine propagation, could then be defined as follows:
 U_{n}^{k+1}=\mathcal{\mathcal{G}}_{\Delta T}\left(U_{n-1}^{k+1}\right)+\mathcal{\mathcal{F}}_{\Delta T}\left(U_{n-1}^{k}\right)-\mathcal{\mathcal{G}}_{\Delta T}\left(U_{n-1}^{k}\right),\quad U_{0}^{k}=u^{0},
 \end{equation}
 with the initial prediction $U^{0}_{n} = \mathcal{G}_{\Delta T}^{n} u^{0}$ for $n=1\ldots N$ and $k=1\ldots K$. $N$ being the number of time subdomains, while $K\geq1$ is the number of predictor-corrector iterations applied. The parareal algorithm is implemented in the library as a separate time integration component, using a fully distributed work scheduling model, as proposed by Aubanel (2010)~\cite{ch5:EA10}. The model is schematically presented in Figure \ref{ch5:fig:FullyDistributedCores}. The parareal component hides all communication and work distribution from the application developer.  It is defined such that a user only has to decide what coarse and fine propagators to use. Setting up the type definitions for parareal time integration using forward Euler for coarse propagation and fourth order Runge-Kutta for fine propagation, could then be defined as follows:

-\lstset{label=ch5:lst:parareal,caption={Assembling a parareal time integrator using forward Euler for coarse propagation and a Runge-Kutta method for fine propagation.}}
-\begin{lstlisting}
+\lstset{caption={Assembling a parareal time integrator using forward Euler for coarse propagation and a Runge-Kutta method for fine propagation.}}
+\begin{lstlisting}[label=ch5:lst:parareal]

 typedef gpulab::integration::forward_euler          coarse;
 typedef gpulab::integration::ERK4                   fine;
 typedef gpulab::integration::parareal<coarse,fine>  integrator;
 typedef gpulab::integration::forward_euler          coarse;
 typedef gpulab::integration::ERK4                   fine;
 typedef gpulab::integration::parareal<coarse,fine>  integrator;

diff --git a/BookGPU/Chapters/chapter7/ch7.tex b/BookGPU/Chapters/chapter7/ch7.tex
index 1e0dd3d7495b1eeffaa76ce728832455e79db10f..a44cbb5aaa8f2411356acd105ff2dc776f57dfa1 100644 (file)
--- a/BookGPU/Chapters/chapter7/ch7.tex
+++ b/BookGPU/Chapters/chapter7/ch7.tex
@@ -470,10 +470,10 @@ The CUDA-based numerical wave model has been developed based on all the numerica
 
 For the unified potential flow model the user will need to provide implementations of the following components; the right hand side operator for the semi-discrete free surface variables \eqref{ch7:FSorigin}, the matrix-vector operator for the discretized $\sigma$-transformed Laplace equation \eqref{ch7:TransformedLaplace}, a smoother for the multigrid relaxation step, and the potential flow solver itself, that reads initial data and advance the solution in time. In order to make the library as generic as possible, all components are template-based, which makes it possible to assemble the PDE solver by combining type definitions in the preamble of the application. An excerpt of the potential flow assembling is given in listing \ref{ch7:lst:solversetup}.
 
 
 For the unified potential flow model the user will need to provide implementations of the following components; the right hand side operator for the semi-discrete free surface variables \eqref{ch7:FSorigin}, the matrix-vector operator for the discretized $\sigma$-transformed Laplace equation \eqref{ch7:TransformedLaplace}, a smoother for the multigrid relaxation step, and the potential flow solver itself, that reads initial data and advance the solution in time. In order to make the library as generic as possible, all components are template-based, which makes it possible to assemble the PDE solver by combining type definitions in the preamble of the application. An excerpt of the potential flow assembling is given in listing \ref{ch7:lst:solversetup}.
 

-\lstset{label=ch7:lst:solversetup,caption={Generic assembling of the potential flow solver for fully nonlinear free surface water waves.}
+\lstset{caption={Generic assembling of the potential flow solver for fully nonlinear free surface water waves.}

 %,basicstyle=\scriptsize
 }
 %,basicstyle=\scriptsize
 }

-\begin{lstlisting}
+\begin{lstlisting}[label=ch7:lst:solversetup]

 // Basics
 typedef double value_type;
 typedef gpulab::grid<value_type> vector_type;
 // Basics
 typedef double value_type;
 typedef gpulab::grid<value_type> vector_type;


@@ -505,10 +505,10 @@ typedef free_surface::potential_flow_solver_3d<potential_flow_types> potential_f
 
 Hereafter, the potential flow solver is aware of all component types that should be used to solve the entire PDE system, and it will be easy for developers to exchange parts at later times. The \texttt{laplace\_sigma\_stencil\_3d} class implements both the matrix-vector and right hand side operator. The flexible-order finite difference kernel for the matrix-free matrix-vector product for the two-dimensional Laplace problem is presented in listing \ref{ch7:lst:fd2d}. Library macros and reusable kernel routines are used throughout the implementations to enhance developer productivity and hide hardware specific details. This kernel can be used both for matrix-vector products for the original system and for the preconditioning.
 
 
 Hereafter, the potential flow solver is aware of all component types that should be used to solve the entire PDE system, and it will be easy for developers to exchange parts at later times. The \texttt{laplace\_sigma\_stencil\_3d} class implements both the matrix-vector and right hand side operator. The flexible-order finite difference kernel for the matrix-free matrix-vector product for the two-dimensional Laplace problem is presented in listing \ref{ch7:lst:fd2d}. Library macros and reusable kernel routines are used throughout the implementations to enhance developer productivity and hide hardware specific details. This kernel can be used both for matrix-vector products for the original system and for the preconditioning.
 

-\lstset{label=ch7:lst:fd2d,caption={CUDA kernel implementation for the two dimensional finite difference approximation to the transformed Laplace equation.}
+\lstset{caption={CUDA kernel implementation for the two dimensional finite difference approximation to the transformed Laplace equation.}

 %,basicstyle=\scriptsize\ttfamily
 }
 %,basicstyle=\scriptsize\ttfamily
 }

-\begin{lstlisting}
+\begin{lstlisting}[label=ch7:lst:fd2d]

 template <typename value_type, typename size_type>
 __global__ void laplace_sigma_transformed(
       value_type* out            , value_type const* p
 template <typename value_type, typename size_type>
 __global__ void laplace_sigma_transformed(
       value_type* out            , value_type const* p


@@ -596,9 +596,9 @@ __global__ void laplace_sigma_transformed(
 
 In a similar template-based approach, the kernel for the right hand side operator of the two dimensional problem is implemented and listed in Listing \ref{ch7:lst:rhs2d}. The kernel computes the right hand side updates for both surface variables, $\eta$ and $\tilde{\phi}$, and applies an embedded penalty forcing \eqref{ch7:eq:penalty}, for all nodes within generation or absorption zones. The penalty forcing functions are computed based on linear or non-linear wave theory in a separate device function.
 
 
 In a similar template-based approach, the kernel for the right hand side operator of the two dimensional problem is implemented and listed in Listing \ref{ch7:lst:rhs2d}. The kernel computes the right hand side updates for both surface variables, $\eta$ and $\tilde{\phi}$, and applies an embedded penalty forcing \eqref{ch7:eq:penalty}, for all nodes within generation or absorption zones. The penalty forcing functions are computed based on linear or non-linear wave theory in a separate device function.
 

-\lstset{label=ch7:lst:rhs2d,caption={CUDA kernel implementation for the 2D right hand side.}%,basicstyle=\scriptsize\ttfamily
+\lstset{caption={CUDA kernel implementation for the 2D right hand side.}%,basicstyle=\scriptsize\ttfamily

 }
 }

-\begin{lstlisting}
+\begin{lstlisting}[label=ch7:lst:rhs2d]

 template <typename value_type, typename size_type>
 __global__ void rhs(value_type const* p    , value_type const* p_surf
                   , value_type const* eta  , value_type* dp_surf_dt
 template <typename value_type, typename size_type>
 __global__ void rhs(value_type const* p    , value_type const* p_surf
                   , value_type const* eta  , value_type* dp_surf_dt

diff --git a/BookGPU/Chapters/chapter8/ch8.tex b/BookGPU/Chapters/chapter8/ch8.tex
index 5d012eb5ec2b5c8047b425b9f43a46de1f0bb72e..f5c001ad7f962bb7e76ed71b663555e9822b4f19 100644 (file)
--- a/BookGPU/Chapters/chapter8/ch8.tex
+++ b/BookGPU/Chapters/chapter8/ch8.tex
@@ -140,7 +140,6 @@ The parallel evaluation of bounds \index{parallel evaluation of bounds} model, a
 \label{ch8:BB-FSP}
 
 \subsection{Definition of the Flowshop Scheduling Problem} 
 \label{ch8:BB-FSP}
 
 \subsection{Definition of the Flowshop Scheduling Problem} 

-\label{ch8:LB-FSP}

 
 As a case study for our GPU-based Branch-and-Bound, we considered the NP-hard and well-known problem in the scheduling theory: the "Permutation Flow-shop Scheduling Problem" (FSP). 
 In this work, the mono-objective case is considered. The FSP aims to find the optimal schedule of n jobs on m machines so that the overall completion time of all jobs, called {\it makespan}, is minimized. 
 
 As a case study for our GPU-based Branch-and-Bound, we considered the NP-hard and well-known problem in the scheduling theory: the "Permutation Flow-shop Scheduling Problem" (FSP). 
 In this work, the mono-objective case is considered. The FSP aims to find the optimal schedule of n jobs on m machines so that the overall completion time of all jobs, called {\it makespan}, is minimized. 


@@ -169,7 +168,6 @@ Figure~\ref{flow-shop} illustrates a solution of a flow-shop problem instance de
 \vspace{0.3cm}
 
 \subsection{Lower Bound \index{Lower Bound} for the Flowshop Scheduling Problem}
 \vspace{0.3cm}
 
 \subsection{Lower Bound \index{Lower Bound} for the Flowshop Scheduling Problem}

-\label{ch8:LB-FSP}

 
 The lower bounding technique provides a lower bound (LB) for each sub-problem generated by the branching operator. The more the bound is accurate, the more it allows to eliminate not promising nodes from the search tree. Therefore, the efficiency of a B\&B algorithm depends strongly on the quality of its lower bound function. In this chapter, we use the lower bound proposed by Lenstra {\it et al.}~\cite{ch8:Lenstra_1978} for FSP, based on the Johnson's algorithm~\cite{ch8:Johnson_1954}.
 
 
 The lower bounding technique provides a lower bound (LB) for each sub-problem generated by the branching operator. The more the bound is accurate, the more it allows to eliminate not promising nodes from the search tree. Therefore, the efficiency of a B\&B algorithm depends strongly on the quality of its lower bound function. In this chapter, we use the lower bound proposed by Lenstra {\it et al.}~\cite{ch8:Lenstra_1978} for FSP, based on the Johnson's algorithm~\cite{ch8:Johnson_1954}.
 


@@ -464,7 +462,6 @@ CUDA enabled devices use several memory spaces, which have different characteris
 The data access optimization challenge is to find the best mapping of the data structures of the application at hand (different sizes and access frequencies) and the GPU hierarchy of memories (different sizes and access latencies). For instance, of these different memory spaces, global memory is the most plentiful but the one with the highest access latency. On the contrary, shared memory is smaller in size but has much higher bandwidth and lower latency than the global memory.
 
 \subsection{Complexity analysis of the memory usage of the Lower Bound }
 The data access optimization challenge is to find the best mapping of the data structures of the application at hand (different sizes and access frequencies) and the GPU hierarchy of memories (different sizes and access latencies). For instance, of these different memory spaces, global memory is the most plentiful but the one with the highest access latency. On the contrary, shared memory is smaller in size but has much higher bandwidth and lower latency than the global memory.
 
 \subsection{Complexity analysis of the memory usage of the Lower Bound }

-\label{ch8:MemComplex}

 
 In this section, the characteristics of the data structures used by the lower bound function are studied in terms of sizes and access frequencies. For an efficient implementation of the LB, six data structures are required: the  matrix $PTM$ of the processing times of the jobs, the matrix of lags $LM$, the Johnson's matrix $JM$, the matrix $RM$ of the earliest starting times of jobs, the matrix $QM$ of their lowest latency times and the matrix $MM$ containing the couples of machines. The complexities of the different data structures are summarized in Table~\ref{ch8:tabMemComplex} where the columns represent respectively the name of the data structure, its size and the number of times it is accessed.
 
 
 In this section, the characteristics of the data structures used by the lower bound function are studied in terms of sizes and access frequencies. For an efficient implementation of the LB, six data structures are required: the  matrix $PTM$ of the processing times of the jobs, the matrix of lags $LM$, the Johnson's matrix $JM$, the matrix $RM$ of the earliest starting times of jobs, the matrix $QM$ of their lowest latency times and the matrix $MM$ containing the couples of machines. The complexities of the different data structures are summarized in Table~\ref{ch8:tabMemComplex} where the columns represent respectively the name of the data structure, its size and the number of times it is accessed.
 


@@ -505,7 +502,6 @@ To reduce the computation time cost of the term $\min\limits_{(i,j)\in \jmath^2,
 \end{table}
 
 \subsection{Data placement pattern of the Lower Bound on GPU}
 \end{table}
 
 \subsection{Data placement pattern of the Lower Bound on GPU}

-\label{ch8:MemComplex}

 
 This section discusses how best to map the six data structures identified above on the various kinds of memories of the GPU device.
 
 
 This section discusses how best to map the six data structures identified above on the various kinds of memories of the GPU device.
 

diff --git a/BookGPU/Makefile b/BookGPU/Makefile
index 22e9fd752bc52d0edc1c2c80b4270aa279744a73..1ab7327c6d5acbb5e1d8725cbe32d237d55026d9 100644 (file)
--- a/BookGPU/Makefile
+++ b/BookGPU/Makefile
@@ -34,5 +34,5 @@ all:
 #      dvipdf ${BOOK}
 
 clean:
 #      dvipdf ${BOOK}
 
 clean:

-       rm -rf *~ *.aux *.bbl *.blg *.dvi *.lon *.not *.lof *.log *.lot *.toc *.ind *.idx *.ilg 
+       rm -rf *~ *.aux *.bbl *.blg *.dvi *.lon *.not *.lof *.log *.lot *.toc *.ind *.idx *.ilg *.aux