/articles/toma

 Using $c_m$, the minimal required bin ($h_m$) is computed as: $h_m \eq floor(h_b + c_m)$, where $h_b$ is the bin accessed by a first thread of the block.  It is computed based on a index of a thread block in the computational grid as described in \tablename~\ref{table:alg_vars}. The $c_m$ is computed during the initialization stage and is stored along with other projection constants in the GPU constant memory. 

 

 

\centering

\includegraphics[width=0.45\textwidth]{img/alutex.pdf}

\caption{\label{fig:alutex} The figure illustrates reconstruction process relaying on the shared memory cache and the algebraic units to perform back-projection. To reconstruct 32x32 square, a thread block caches 48 bins from each projection row. The projections are processed in groups moderated by the size of available shared memory. First, the required subset of bins in each projection is determined (left). The selected subsets along with their offsets in the projection rows are cached in the shared memory (center). Then, the reconstruction is performed and projections are processed in a loop one after another (right). Each thread is responsible for several pixels of output slice. For each pixel the required position in the sinogram is computed. The cache offset for the considered projection row is subtracted from this position and the offset in the cache is determined (bottom). As the offset is typically not integer, two array elements are loaded from the cache and interpolation is performed.}

 

  In the first stage of algorithm, the number of threads assigned to cache each projection is adjusted to optimize access to the shared memory. If a large 64-by-64 area is reconstructed, a full warp of 32 threads can be assigned for each projection row avoiding any possible bank conflicts. Unfortunately, it is not completely optimal on the Kepler architecture as, then, it is impossible to re-combine two bins into a single 64-bit shared memory write as explained above. It is also not possible to assign 32 threads per row for a smaller 32-by-32 area because only 48 bins has to be cached per projection in this case. And it is a bad idea to keep the half of the threads idling part of the time. Therefore, several projection rows are processed by each warp in the described cases. This potentially may result in the bank conflicts. If only a single slice is reconstructed, however, the banks are shifted from one projection row to another as illustrated on \figurename~\ref{fig:banks}. On the platforms with 32-bit shared memory, the caching is performed without causing bank conflicts if groups of 16 threads are used per projection row. On the Kepler platform, however, only 8 threads per projection row are used instead and each thread is assigned to write two values into the cache. 16 threads are used to process 32-bins per iteration on the Kepler if 64-by-64 area is assigned to thread block. For the multi-slice reconstruction modes, 16 threads per projection are optimal on all platforms. On Kepler, the \emph{float2}-sinogram used together with the 64-bit banks is performing exactly the same as the \emph{float}-sinogram is performing with the 32-bit banks. I.e. 16 threads per projection are required to avoid bank conflicts. On other platforms it is enough to prevent bank conflicts within the groups of 16 threads while dealing with 64-/128-bit data. So, there are no problems if 16 threads are used per projection row. Two shared memory buffers, however, are used as explained above if 4-slice reconstruction is executed on  any AMD device or Fermi-based GPU from NVIDIA. The optimal settings for each reconstruction mode are summarized in \tablename~\ref{tbl:shmemconf}.

  

\centering

\includegraphics[width=0.45\textwidth]{img/banks.pdf}

\caption{\label{fig:banks}  The figure illustrates how the warps are assigned to cache a subset of a sinogram on the systems with 32-bit and 64-bit shared memory. For each projection 48 bins required to reconstruct area of 32-by-32 pixels are cached. The shared memory banks used to back each group of 16 bins are specified considering that 32-bit data format is used.}

 

  In the reconstruction stage, a number of bins accessed by each warp should be minimized. If the linear mapping is used, up to 32 projection bins may be accessed by a warp. This has no negative impact on the performance in the single-slice mode. On Fermi architecture, however, the shared memory bank conflicts are possible if 4 slices are reconstructed in parallel. To prevent conflicts, it is necessary to restrict the loads by the half-warp to 8 bins only (which spread over exactly 32 banks if \emph{float4} data is used). This can be ensured if half-warps are mapped to the pixel squares of 4-by-4 pixels each. The  256 threads of the block are mapped to 16 such squares. For the reasons explained in \sectionname~\ref{section:newtex}, the squares are arranged linearly along $x$-axis. Two rows of 4x4 squares are processed in parallel if a small 32-by-32 area is reconstructed. Only a single row is covered for the bigger area or if 128 threads per block are used. The remaining rows are processed over 4/16 iterations as it was done in the texture-based kernel. The threads accumulate the sums for each pixel in the register-bound array and dump it to the global memory once processing of all projections is completed.

 

\centering

\includegraphics[width=0.45\textwidth]{img/alumap.pdf}

\caption{\label{fig:alumap} The assignment of block threads to pixels as proposed for ALU-based reconstruction }

\input{table_52_shmem.tex}

 

\subsubsection{Advanced Caching Mode}\label{section:ycache}

  

\centering

\includegraphics[width=0.45\textwidth]{img/ycache.pdf}

\caption{\label{fig:ycache} Advanced Caching Mode. Two values are cached for each bin of a sinogram. The second value stores the difference between neighboring bins to allow faster interpolation. The shuffle instruction is used to get values from the bins cached by a different GPU thread.}

 To cache bins required to reconstruct 32-by-32 pixel area, we use 16 threads per projection row on both 32- and 64-bit platforms. The banks are shifted between projection rows on the 64-bit platform as explained in \sectionname~\ref{section:alurec_shmem}. And on 32-bit system it is enough to avoid bank conflicts within a half-warp. Furthermore, there are also no bank conflicts between the threads of a half-warp as the stride is not multiple of 4. See illustration in \figurename~\ref{fig:ybanks}. A full warp is used per row if a thread block is assigned to process larger 64-by-64 pixel area. The same number of iteration is, then, required to process the complete projection row and, consequently, the shared memory is accessed with the same stride without bank conflicts.

 

 

\centering

\includegraphics[width=0.45\textwidth]{img/ybanks.pdf}

\caption{\label{fig:ybanks} The figure illustrates how the shared memory banks are accessed if advanced caching mode is used, see \sectionname~\ref{section:ycache}. The presented layout is employed to reconstruct area of 32x32 pixels. The grayed boxes indicate the banks accessed by a warp during the first caching iteration. Two values are cached for each bin and, consequently, each bin spans over two memory banks on the platforms with 32-bit shared memory. On these platform it is also only necessary to avoid conflicts within a half-warp. So, the accesses for second half-warp are not shown. }

 

 We evaluate a single-slice reconstruction ($n_v \eq 1$) using advanced caching mode. The blocks of 256 threads ($n_t \eq 256$) are assigned to process 32-by-32 pixel square ($n_q \eq 4$). For the sake of simplicity we assume 16 threads are used to cache a single projection row and that the inner projection loop is fully unrolled. Then, the following number of operations is required per single update. We skip the texture fetches as load is very low and certainly is not a limiting factor. To verify our estimates, we run the prototype implementation under CUDA profiler. The number of estimated and measured operations is compared in \tablename~\ref{tbl:modelprof}. There is difference, but the error is within 10\%.

 

\begin{threeparttable}

\caption{\label{tbl:modelprof} Estimated and measured number of different operations required to perform back-projection using linear interpolation}

\centering

\item \textbf{Instructions}: 8.2

\end{itemize}

 

 

 As can be seen the performance bottleneck is architecture dependent. The AMD VLIW and NVIDIA Kepler GPUs are bound by ability to perform rounding operations and convert variables between integer and floatint-point representation. While not limiting performance in the measured configuration, it still sets quite low threshold on Maxwell and Pascal GPUs. However, the main limiting factor on these architectures is shared memory bandwidth. The Fermi GPU is only capable to dispatch a single instruction per warp and, consequently, bound by the ability instruction throughput. Finally, AMD GCN based devices are restricted by the performance of algebraic units.

 

\begin{threeparttable}

\caption{\label{tbl:model} Performance estimates according to performance}

\centering

\end{lstlisting}

 

 

\centering

\includegraphics[width=0.45\textwidth]{img/fancy.pdf}

\caption{\label{fig:fancy} IEEE 754 representation of single-precision floating-point number (top) and an example how to get the standard integer representation in fraction part by adding $2^{23}$ (bottom)}

 

The performance along with number of issued instructions per update is shown in \tablename~\ref{tbl:fancy} for NVIDIA GTX Titan. The speed-up of 20\% is achieved if rounding is implemented using floating point instructions, but index computation is left on SFU. The complete elimination of SFU instructions puts unnecessary load on ALUs keeping SFU idling. The difference between the two optimization modes is hard to see in profiling results since integer multiplications are just replaced with integer additions. The last operation is performed by ALU, but both are counted as integer operations by NVIDIA profiler. This method has little impact on VLIW architecture. While the performance is limited by the throughput of integer instructions, the difference between performance of floating-point and special units is not as high as on Kepler. Consequently, performance is limited at  the same level approximately.

 

\begin{threeparttable}

\caption{\label{tbl:fancy} Different interpolation modes on Titan (Kepler) GPU}

\centering

 \item The number of threads assigned to cache each projection row is in accordance with the requirements specified in \tablename~\ref{tbl:shmemconf}. Then, no shared memory bank conflicts happen in the caching stage of reconstruction and the shared memory performance writes are executed optimally.

 \item The number of projections cached at each iteration step is enough to utilize all threads in the block. Otherwise, the threads idling at the synchronization point reduce the efficiently achieved occupation.

 \item The projection loop is completely unrolled to provide additional ILP parallelism and ensure that multiple 32-bit memory operations can be combined into a single 64-/128-bit instruction.

 \item All appropriate optimizations discussed through this section are implemented.

\end{itemize}

 

 

 The Maxwell and Pascal GPUs have a large amount of shared memory and registers, but are bound by the shared memory bandwidth. Since the relative amount of the shared memory writes is reduced for the larger pixel area, the 64-by-64 pixels are processed by each thread block on these platforms. In the linear interpolation mode the amount of the shared memory operations is well balanced with ALU throughput.  The streaming of memory reads is not required if the shared memory loads are interleaved with ALU- and SFU-bound interpolation instructions. Thus, the 100\% occupancy can be reached. The speed-up of 15\% is measured compared to the default compiler settings if the desired occupancy is not hinted.  This is not the case in the nearest neighbour interpolation mode. The shared memory bandwidth is the bottleneck in this case and the performance is improved if multiple shared memory loads are streamed together. Consequently, significantly more registers are required. By default the CUDA compiler does not utilize the streaming capabilities fully, but runs at 62\% occupancy. Requesting even lower occupancy of 50\% allows to stream more loads together and improves performance by 7\%. The impact of occupancy on the performance for both linear and nearest-neighbour interpolation modes are reviewed \tablename~\ref{tbl:occupancy}.

 

\begin{threeparttable}

\caption{\label{tbl:occupancy} The effect of occupancy-targeting on the performance for NVIDIA Titan X GPU}

\centering