\subsection{What Makes Convolution Parallel Friendly?}
Before attempting to parallelize an algorithm it is good to know that it is parallelizable. Convolution is an imaging technique, and like most imaging techniques one is performing similar actions on several pixels \cite{Gonzales}. If these actions are independent as in the case of convolution, then one can run actions of each pixel in parallel with other neighboring pixels. This is actually part of the motivation behind the structure of the GPU and the way it is used (considering its SIMD architecture). Convolution also meets all the requirements of Bernstein's conditions (as mentioned in Section \ref{parallelsection}) in that each pixel is completely independent of the others. Another interesting aspect in convolution is that calculations  performed per pixel are actually quite time consuming in that each pixel must travers the length of the filter mask. The fact that being able to mask other computations behind each calculation should show great speedup.  These are some of the motivation behind parallelizing convolution and why it is such a parallel friendly technique. The major bottle neck here is that since so many calculations are performed per pixel is the correct use of caches such that calculations are done in an optimal fashion.

\subsection{Spatial vs. Fourier Domain}
This project has a scope only accounting for spatial domain filtering as it is the most expensive of the two and because of time limitations. Convolution can be perormed in the fourier domain, but has its limitations. For example, The fourier domain is asyptotically faster than the spatial domain, but has limitations such as the data set dimensions have to be a power of two \cite{victor}. This is usually not the case for seismic data, and therefore one would have to pad the data, which could cause distorsion or noise in the filtered image. This gives justification that performing convolution in the spatial domain is also usedfull in some cases, and to account for the limitations of the time available, our project focuses on the spatial domain.

\subsection{Memory Access Patterns in Uni-Core CPU Implementation}
In the first implementation only the use of one CPU core was in focus. This of course has the advantage of having to narrow the focus on optimizing memory use for one core, which then should be easy to expand to the use of quad cores with the use of threads, more precisely in this case the use of OpenMP. 

Lets have a closer look at how the memory works in this case. After having read to the buffer there are three levels of cache that can be utilized \cite{amd}. The data in the buffer is aligned according to rows read from disk. This is a problem because when filtering for a pixel we will be looking at all its neighbors in all three dimensions. Which means that one must perform several jumps in memory to obtain the next values. This is of course not ideal, but since the jumps in memory are regulated by the dimensions of the filter, which are constant, the pre-fetcher should be able to recover them while calculations are being perform resulting in close to no cache misses. This is seen in the profiler results where when running on one CPU the cache misses are 0.01\% of the time. Because of the nature of how the data is aligned and the constant jumps in memory the prefetcher is actually able to memory optimize the convolution algorithm. 

This is the case of fetching from the buffer to the third level of cache. In the case of when the correct data is retrieved in the third level of cache. One is then retrieving a subset of that data to the level 2 and level 1 cache for calculations. The nice thing here is that now the data is aligned such that each sub cube, of an equal size as the filter, are placed after one another. And when performing calculations in the level 1 cache then the data needed is aligned such that it is read in a straight line resulting in well aligned data on all levels of  the architecture. The only cache misses present are when one reaches the end of a the data set in the x dimension, and proceeds to jump to the next row below it. This cache miss occurs on the level 2 and level 1 cache and is a result of non-overlapping data between the two pixels being filtered. In the normal case there is overlap between neighboring pixels and one can perform calculations on these values since they are in the cache while the pre-fetchers attains the next values, but in the case where there is no overlap one must wait for the values to be gathered before proceeding. This is a rare case, because of the blocking pattern, and that is why it does not effect the overall cache misses of the application.

\subsection{Memory Access Patterns in Quad-Core CPU Implementation}
Knowing that one has an optimal serial implementation one can then justify comparisons with other implementations to analyze gained speedup. In the quad-core implementation one can use this optimal structure of memory and aligned data to run 4 threads where all the data in the level 3 cache is shared. This should result in that the implementation be more bound by memory bandwidth since memory access is optimized. Given that there are now 4 CPU performing the same task as one did before and are totally independent this should result in a theoretical 4 time speedup. Some bottlenecks here are of course the memory transfer rates, but the implementation should come quite close for large amounts of data and computation. In this project, the results obtained are a 3.57 speedup, which is close to the theoretical value and also underlines the fact that the implementations are memory optimized. The profiler also indicates this in the same way as for the uni-core implementation. 

\subsection{Effects of Operating System on Convolution}

In the case of running on the uni-core one has 3 other cores that the operating system can use to perform tasks other than that of convolution. But, in the case of utilizing all the capacity of the CPU by running the convolution process on all 4 cores, tasks related to the OS will interrupt or interleave with the convolution process. This could be one of the major issues related to the speedup limitations noticed in the results. The fact that a speedup of 3.57 and not an ideal 4 times speedup is achieved when running on several cores could have many reasons. One could be the fact that  only convolution is run in parallel and not I/O. But the trend that even for large sets, where I/O is less significant, no greater speedup was seen could be because large sets have greater execution time and the operating system is bound to interrupt the process, creating another bottleneck. This is avoidable by running on an accelerator such as the GPU. 

\subsection{Effects of the PCI-express Bus}
Even though the same operation is performed on the GPU the memory access patterns differ. First of all there is a an addition to the pipeline, which is actually quite time consuming \cite{nvidia_cuda_prog_guide}, and is the transfer of data from the system buffer to the GPUs buffer (global memory). This is done on the PCI express bus which has an average bandwidth of about 1 GB per second. This can prove to be troublesome if one has to perform several transfers such as in the case of large data sets. Another issue is that on the GPU there will be more transfers because of the size of the memory. The system buffer has a size of 8 GB ,where 2 are occupied by the OS. This allows 3GB (since convolution requires duplication) to be read from disk at a time as discussed earlier. While the GPU memory only allows 2GB. Therefore the combination of both having to transfer more often to process the same amount of data and that the transfer is slower because of the bottle neck the PCI express bus one cannot expect linear speedup between the amount of cores on the GPU and the CPU. Not to mention that the GPU processing units have a lower frequency. Nevertheless given that convolution is very parallelizable a speedup is expected. 

\subsection{Results of CPU vs GPU}
The results indicate that for large data sets that the GPU is actually about 17 times faster than a uni-core CPU. And is expected to be faster as long as each thread has a lot to do i.e. the size of filter is large. Because each pixel has to calculate values along the filters dimensions. Nonetheless a 17 time speedup is actually a good result in the case of pre-processing large amounts of data. In the case where one would use 16 hours on a uni-core. On the GPU it would run in about 54 minutes. Thereby one has reduced the problem dimensions from hours to minutes.

\subsection{GPU Occupancy}
To utilize the processing capabilities of the GPU one has to make sure that there are enough threads available to actually perform all the calculations and obtain maximum occupancy on the GPU. NVIDIA have produced an occupancy calculator \cite{nvidia_cuda_occ} that helps developers to analyze their implementations and show how much of the GPU is actually occupied  at a time. Here there is a fine balance between the number of threads, blocks and registers one uses. Because there is a limited number of registers per kernel then restrictions are introduced to the max amount of threads. In the case of this convolution algorithm there is a need for 23 registers (information retrieved from the cubin file), which results in only 63\% of possible threads can be run at a single kernel launch. Meaning that only 63\% occupancy can be obtained. The alternative is to push some values from the registers to global memory and run more threads. But since the register values are used often and forced lowering of register use in the complier would be place values in global memory will result in lower performance. This is because registers are accessed much faster than global memory and the values in these registers are used constantly.

\subsection{Memory Access Patterns in GPU Implemenatation}
Other than transferring the data to the graphical processor one has to optimize the kernel (algorithm running on the GPU) such that it utilizes the performance capabilities of the GPU. As mentioned before, convolution is a very bandwidth bound problem and on the GPU there are several layers of memory that can accelerate bandwidth, but they need to be explicitly programmed to. The shared memory layer is a private memory for each block such that threads within the block can communicate and cooperate without having to access global memory. This is memory that can be accessed within one cycle and is on chip, but needs to be transfered from global memory \cite{nvidia_cuda_prog_guide}. In the case of convolution there is no communication or cooperation within the threads and therefore it is not used. There is another layer of memory that can prove to be useful in this case, and that is constant memory. The issue here is that constant memory is small in size. One could use the global memory to run ones code from, but the memory latency would ruin performance \cite{nvidia_cuda_prog_guide}. 

\subsection{Why use Constant Memory and Not Shared Memory?}
Shared memory is the fastest of all buffer layers on the GPU. It can be accessed within one cycle of computation. The only other buffer that has this attribute is constant memory, and even here it is only valid for best case when the data is coalesced. The trouble with both these buffer is their size. Both are in the range of KB, even though constant memory is 4 times larger. The issue with shared memory is that it needs to be updated per kernel launch and one cannot send data directly to it since it is on chip. To use shared memory one has to retrieve data from the global memory and then do calculations on them. The trouble is that since it is so small (16KB) not much data will be buffered there and the overhead work of transferring the data from global memory might result in a decrease in performance. Making its use both complicated and with no guarantee of speedup.

Constant memory on the other hand is off-chip and can be updated directly without having to go through global memory, and it does not require to be updated per kernel launch. But, it must be allocated statically. Meaning that one must at compilation time decide how much memory should be allocated. This is ideal for values that are constant throughout the computations, hence its name. Another limitation here is that constant memory is read-only from device. Therefore it cannot be used to store convoluted data, but in this implementation it is used to store the convolution filter. This is ideal since the convolution filter is used in all the calculation and retrieving results from it is half the job in convolution. Therefore if coalesced it can be retrieved in one cycle (best case) resulting in great speedup. And since the filter values are read only this is ideal. If the filter is to be placed in shared memory then it would have to retrieved from global per kernel launch, which is a great overhead compared to the use of constant memory. Shared memory also has to be coalesced to perform well.

\subsection{What if Convolution Filter is Larger Than Constant Memory?}
For filter sizes larger than 25x25x25 there will not be enough space in constant memory. There are two solutions here either one can run everything from global memory, which the results show gives a 12 times speedup. Or, one could perform several iterations with a smaller filter. The result would be visually similar, in the case of blurring, and computationally there are similar amounts of operations performed. Mathematically however it is not similar and a larger filter should be used. A larger filter would only result in more computations, but the trouble with using smaller filters is that one would double communication time and time spent transferring data to the GPU, since the operation is done several times. Therefore placing all data in global would most likely give a better execution time. Since disk access and memory bandwidth are the main bottlenecks and doubling their time of execution would probably result in a sub-optimal solution.

\subsection{Bank Conflicts and Coalescing}
Another aspect of memory usage on the GPU that is worth noting is the possibility of bank conflicts that happen during writes to the memory. The way this is solved here is that each thread solves for a value and that is why each neighboring pixel writes in the neighboring address in memory. This results in coalesced writes, which are recommended by NVIDIA \cite{nvidia_cuda_prog_guide} to optimize performance. Not only does the GPU algorithm have coalesced writes, but it also reads coalesced. This is because if the current thread is reading a value, its neighboring thread is reading the value beside it at all times. Again this is recommended for high performance. 

\subsection{The Use of Multi-GPU}
Now that all the steps are taken such that the kernel is running optimally on a single GPU one could progress to running the code on several GPUs. The results show that running on 4 GPUs gives a speedup that is 3.7 times faster than on a single GPU for larger sets of data. This is the result of running the kernels from 4 threads on the host where each invoke their own GPU, and the pipeline structure in which data is transfered to run is also optimal. The fact is given the way that the application executes on the multiple GPUs, with running GPUs one after the other and not interleaving at the data transfer stage, hides a lot of computation and transfer time. And as the problem size grows the more evident this becomes. Because the only delay evident that will not be there when running on a single GPU is the three cudaMemCpy instructions run at the end of the execution. For large data sets this takes some time, but is not a very heavy task that is why one can obtain a speedup of close to 4. This is even more evident in the reduction of computation time as the data size grows for the 4 GPUS. since when there is little data, calculations will be fast and the transfer of data to and from the GPU will be the bottleneck. While in the case where there are larger amounts of data, the calculations are larger and the time of the transfer of data is smaller in comparison to calculations. The transfer of data to and from the GPU is counted as computation time since it is an overhead of using the GPU contra the CPU, and since its a new introduction to the computation pipeline then it should be accounted for in computation. 

\subsection{Communication vs. Computation time}
One of the largest motivations to parallelize the convolution algorithm is the amount of time that is spent on computations rather than retrieving the data from disk. When running on a single CPU 99\% of the time was used for computations. When trying to run the algorithm on several cores the amount of computation time was reduced to 94\% on the quad-core and 75\% on a single GPU. It was not before running on multiple GPUs a 10\% computation time was achieved. This shows that convolution as a problem scales well with several cores. This is because of the independence of calculations. When reaching a distribution of 10\% computation time and 90\% I/O time then one should realize that further paralleizing and optimizing the 10\% computation time will result in little gain, and the focus should now be shifted at improving I/O access. Maybe the use of parallel I/O could be useful, but given the architecture of the machines used here there is only one node (CPU) able to access disk and buffer. This means that there are no options to run disk access in parallel, but it would be an interesting aspect for further work.