As mentioned before, the GPU is an accelerator and therefore it must cooperate with the CPU. This introduces some delays in sending the data from disk to the global memory of the GPU, Since it has to go through the main memory and be streamed over to the GPUs memory. There is also a delay after the computation is done since the data has to be sent back to main memory from the global memory of the device (GPU). This means that there is some extra work that has to be done in order to perform computations and may result in no speedup in the case where little computation is needed. The communication between CPU and GPU is done across a PCI-express buss which has a limited bandwidth of about 6GB/s, which is less than optimal. But, since we are looking at large data sets and a computationally heavy algorithm this communication time can be shadowed by computation time. Meaning that the use of the GPU should be good for performance in this case.

When using the GPU for calculations one uses several threads running on the device to perform calculations. There is a time constraint of 5 seconds that should not be exceeded by a kernel invocation, which means that the tasks performed per thread should be of a simple nature. In this case each thread can be performing calculations on one pixel of the data cube at a time. This should be a simple task even for a larger filter. Here similarly to when running threads on the multi-core CPU the memory is shared by all the threads. Thereby making it easy to allocate global memory, transfer the data there, invoke the kernel and get the results. 

The main bottleneck in this procedure is the transfer of data to and from the GPU since it is limited by the capacity of the bus and is surely slower than the rate data can be retrieved from the global memory for calculations on the GPU, 102 GB/s \cite{nvidia_cuda}.  That is why it is a good idea to transfer as much data as possible at a time and then invoke the kernel several times to perform all the necessary calculations. Here as well we have the issue that convolution needs a equally large array to store its results therefore forcing the use of only half the memory for data transfer at a time. When it comes to the memory usage internally on the GPU one has several options. One could allocate the data in the shared memory, which is useful if one is to reuse data between threads. Or one can use the texture memory, which is optimized for two dimensional texture calculations. In this implementation none of these intermediate memory banks are used because of the lack of communication between threads given the nature of the task. The threads do not need to share memory in order to solve their tasks since there are no dependencies. The filter, on the other hand is located in constant memory since it will not be changed during the process and such that all the threads can read it. It is smart to use the constant memory becaused it is actually cached. The only limitation is that it is read only by the device. This is not an issue because one is only interested in reading the filter values and not change them during execution. Another interesting aspect that should be noted is the fact that the best case lookup in constant memory is one cycle, which is similar to shared memory.  

\begin{figure}[H]
	\centering
	\includegraphics[width=12cm]{report/method/gpucoal.PNG}
	\caption{Image illustrating memory coalesced reads \cite{nvidia_cuda_prog_guide}}
	\label{coalesced}
\end{figure}

In this implementation the amount of threads created per block has been decided with the use of the occupancy calculator developed by NVIDIA \cite{nvidia_cuda_occ}. To use the calculator a cubin file was compiled with the use of the CUDA file containing the convolution kernel. The cubin file can be opened as a text file and contains relevant information such as the use of shared memory and the amount of registers used. After inserting this information into the occupancy calculator it indicated an amount of 100 threads would give the most occupancy, and was at 63\%. The bottle neck here is the amount of registers available per thread, and since one needs 24 to run the kernel one cannot run it on all the cores available on the GPU. The reason being the maximum amount of registers available in total. One could force the compiler to use only 16 threads, but this will force the threads to save intermediate data to global memory for calculations, which will greatly decrease performance. Given that 100 threads would give the maximum occupancy in this case means that one can have several blocks. Since the data is in three dimension the grids two dimensional nature was combined with one dimensional threads to access all three dimensions in the problem. This is a similar approach to the one performed in \cite{aksnes}. 

To copy data over and from the global memory of the GPU standard CUDA functions where called, and even when allocating memory on the main memory on the host CUDA non device functions where used. The reason for this is because it optimizes the transfer and locality of the data read from the disk to match that of the GPU device. This is because of the two dimensional nature of the device memory \cite{nvidia_cuda_prog_guide}. When the data is in global memory and the threads are calculating each their pixel the data locality, by default, will have them avoid writing to the same address in memory avoiding bank conflicts and it will also have them read data coalesced. Since each neighboring thread is working on the neighboring pixel in the result array they will never write to the same value and avoiding bank conflicts. And since each thread will start to read values with an offset of 1 from its neighboring thread they will be reading data coalesced.

Looking at the big picture, the data is now read from the disk to the main memory and then moved to global memory on the device before performing calculations on the GPU. Then the convolution kernel is called as many times as needed given the size of the sub cube in the x dimension (this is because using 100 threads is optimal \cite{nvidia_cuda_occ}). Then data is written back to main memory and further on back to disk before the next block of data is read. The GPU used in this implementation has 4 GB of memory and so only 2GB of memory are read at a time because of the result table. This means that there are several reads from disk needed here compared to the CPU version where one could read 3GB at a time. This could be a possible source of delay since reading from disk and transferring data to memory is time consuming. 
\newline

\lstset{frame=single, tabsize=2, language=Java, basicstyle=\footnotesize, numberstyle=\footnotesize, numbers=left, captionpos=b, breaklines=true, title={Kernel code for three dimentional convolution}}
\begin{lstlisting}
__global__ void increment_kernel(float *d_data, float *d_token, int len, int x, int y, int offsetnrx, int offsetnry) { 
	float imagevalue =0;
	for (int l = 0; l < len; l++) {
		for (int m = 0; m < len; m++) {
			for (int n = 0; n < len; n++) {
				float temp = d_token[(blockIdx.y+l)*( x + len)*(y + len) + (blockIdx.x+m + (gridDim.x*offsetnry))*( x + len)+ (threadIdx.x+n +(blockDim.x*offsetnrx))];
				float temp2 = d_filter[(len*len)*l + (len)*m + n];
				imagevalue += temp * temp2; 
			}
		}
	}
	d_data[(blockIdx.y)* x * y + (blockIdx.x +(gridDim.x*offsetnry))* x + (threadIdx.x)+(blockDim.x*offsetnrx)] = imagevalue / (len*len*len);
}
\end{lstlisting}

