In the case of using multi GPUs one is trying to hide computation time by running computation simultaneously on all GPUs. This is of course not quite possible because of the restrictions of the PCI bus connection when communicating and transfering data to the GPUs. In this implementation the system running has 4 GPUs available with 4 GB of memory each. The code to run and transfer data to the GPUs is run in threads with the help of OpenMP just like in the multi core CPU version. Here similarly to the quad-core implementation there are 4 threads running where each thread is responsible for a GPU and they have corresponding IDs. 

The idea is to run each thread in parallel to invoke all the GPUs simultaneously. The only issue here is how the threads should share the bus. This needs some clever pipelining. In this implementation one thread at a time transfers data (both to and from the device) to utilize the buses capacity and not interleave data on the bus. Then by running the kernel invokation in parallel one can hide the transfer times for the 3 other GPUs since the first GPU is already doing calculations while the others are transferring data. This way the only overhead compared to running on a single GPU would be the the final transfers from the last 3 GPUs, which occur after the first GPU is done computing and transferring the data back. Here all reads and writes are done by a single thread and are as large as possible. To make sure that each thread was accessing the bus separately the use of the ''critical section'' directive in OpenMP was used. This works as a semaphore and avoids deadlocks and collisions in memory access. See Figure \ref{multigpupipeline} for details around the pipeline. It is worth mentioning that the larger the calculations the less signinficant the cudaMemCpy() functions will be in terms of delay. Therefore the size of each block in the figure is not representative for larger sets, but should give an idea of the concept of hiding calculations by runnig the code in threads.

\begin{figure}[H]
	\centering
	\includegraphics[width=14cm]{report/method/multigpu.PNG}
	\caption{Pipeline representing running code on 4 GPUs}
	\label{multigpupipeline}
\end{figure}

On the machine running the 4 GPU solution there is only 12 GB of main memory available, which should be an issue since the quad GPUs together have 16 GB of global memory. In theory one could not transfer all the data from one read of disk to memory to device memory. But, in the case of convolution and the need for a result table, one only needs to transfer 8 GB of data at a time and the rest should be allocated for the result array. Therefore in this case one read from disk would suffice to utilize all the memory available on the graphical processors. 

Another issue worth mentioning is that there is a restriction in CUDA that if several devices are invoked. A CUDA \textit{non-device} function cannot be called between every time one sets the device i.e chooses a device.  The use of non-device functions, such as cudaMallocHost(), help because they optimize reads from host to device and not being able to use them means that one should use standard methods of allocating memory. This is a potential cause for less speedup because the memory is not structured for optimal transfer, but this is just one of the limitations of CUDA. One could potentially try to restructure the data without using the functions to utilize for the GPU.
\newline

\lstset{frame=single,tabsize=2, language=Java, basicstyle=\footnotesize, numberstyle=\footnotesize, numbers=left, captionpos=b, breaklines=true, title={C code with OpenMP for three dimentional convolution on Multiple GPUs}}
\begin{lstlisting}
#pragma omp parallel default(shared)
{
		float *d_data;
		float *d_token;
		int th_id;
		int newtoken;
		int newdata;
		dim3 threads = dim3(100, 1, 1);
		dim3 blocks  = dim3(100, z/num_th, 1);
		
		th_id=omp_get_thread_num();
		cudaSetDevice(th_id);
		int device;
		cudaGetDevice(&device);

#pragma omp barrier

		cudaMalloc((void**)&d_data, d_dataLength);
		cudaMalloc((void**)&d_token, d_tokenLength);
	
		newtoken = (d_tokenLength/4)*th_id;
		float *temptoken = (float*)(token + newtoken);
	
		cudaMemcpy (d_token, temptoken , d_tokenLength, cudaMemcpyHostToDevice);
		increment_kernel(d_data, d_token, len, x, y, blocks, threads);

#pragma omp critical
{		
		cudaMemcpy (tempdata , d_data, d_dataLength, cudaMemcpyDeviceToHost);
}		
		cudaFree(d_data);
		cudaFree(d_token);
		cudaThreadExit();
	}
}
\end{lstlisting}