Concurrent programming with CUDA. Part 2: GPU hardware and parallel communication patterns

Content

Parallel Communication Templates

• run on a stream for each output element, where each stream will average the values ​​of 3 neighboring input elements - gather ;
• or run a stream on each input element, where each stream will add 1/3 of the value of its input element to the value of the corresponding output element - scatter .

1. All the data needed for one stream is grouped in memory (in the case of a one-dimensional array, it is a solid “piece” of memory, in the 2D case, several pieces of memory located at the same distance from each other).
2. The value of some input element is read several times from neighboring streams (the specific number of readings depends on the selected mask) - it becomes possible to "reuse" the data provided to us by CUDA - this will be discussed later in the article.

GPU hardware

• A CUDA-compatible GPU consists of several (usually dozens) streaming multiprocessors (streaming multiprocessors), hereinafter referred to as SM .
• Each SM , in turn, consists of several dozen simple / streaming processors (SP) (regular / stream processors), or more precisely, CUDA cores ( CUDA cores ). These guys are more like the usual CPU - they have their registers, cache, etc. Each SM also has its own shared memory (shared memory) - a kind of additional cache that is accessible to all SPs, and can be used both as a cache for frequently used data and for “communication” between threads of one CUDA block.
• The GPU also has its own memory, called device memory , common to all CUDA streams - it is with it that the cudaMalloc and cudaMemcpy functions work (also schoolchildren and GPU manufacturers like to measure it by size) .

Compliance with the CUDA model and GPU hardware. CUDA Warranties

• Each block will be fully executed on the allocated SM .
• The distribution of blocks by SM is a GPU, not a programmer.
• All flows of block X will be divided into groups called warps (usually they say warps) and executed on SM . The size of these groups depends on the GPU model, for example, for models with Fermi microarchitecture it is 32. All threads from one warp are executed simultaneously, occupying a certain part of SM resources . Moreover, they either execute the same instruction (but on different data), or are idle.

• All threads in a specific block will be executed on any one SM .
• All threads of a particular kernel will be executed before the next kernel starts.

• Some block X will be executed before / after / simultaneously with some block Y .
• Some block X will be formed at some specific SM Z .

Synchronization

• A barrier is a point in the kernel code, upon reaching which a thread can "go" further only if all threads from its block have reached this point. Once again: the barrier allows you to synchronize only the flows of one block , and not all flows in principle! The restriction is quite natural, because the number of blocks set by the programmer can significantly exceed the number of available SMs .
• Atomic operations - similar to atomic operations of the CPU, a complete list of available operations can be found here .
• __threadfence is not quite a primitive of synchronization: upon reaching this instruction, a thread can continue execution only after all its memory manipulations become visible to other threads - in fact, it forces the thread to execute flush cache.

Key Principles for Using CUDA Effectively

• The principle of increasing the ratio (time of useful work) / (time of operations with memory) was discussed in a previous article. The value of a fraction can be increased in two ways - increase the numerator, decrease the denominator: that is, you need to either do more work or spend less time on memory operations. In addition to the obvious solution - to reduce the number of memory accesses as much as possible, the following principles of effective work with memory are used:
  
  • Moving frequently accessed data in the faster memory: thread local storage > total memory block >> shared memory device >> memory of the host. Thus, if several threads use the same data in the same block, it most likely makes sense to move them to the shared memory of the block.
  • Sequential access to memory: since the flows in blocks are actually executed in warp groups, provided that the threads in one warp work with data located in memory sequentially, CUDA can read one large piece of memory in one instruction. Otherwise, if the streams in the warp access data scattered in memory, the number of memory accesses increases.
  
  
• Reducing the divergence of flows: a feature of CUDA is that flows in one warp always either execute the same instruction or are idle. Thus, if there is a code in the kernel code of the form
  

  if (threadIdx.x % 2 == 0)
  {
    ...
  }
  ...
  

  
  , then half of the warp threads (those with an odd index) will wait until the second half executes the code inside the if-a. Therefore, such situations should be avoided.

We write the second program on CUDA

#include<chrono>#include<iostream>#include<cstring>#include<string>#include<opencv2/core/core.hpp>#include<opencv2/highgui/highgui.hpp>#include<opencv2/opencv.hpp>#include<vector_types.h>#include"openMP.hpp"#include"CUDA_wrappers.hpp"#include"common/image_helpers.hpp"voidprepareFilter(float **filter, int *filterWidth, float *filterSigma){
  staticconstint blurFilterWidth = 9;
  staticconstfloat blurFilterSigma = 2.;
  *filter = newfloat[blurFilterWidth * blurFilterWidth];
  *filterWidth = blurFilterWidth;
  *filterSigma = blurFilterSigma;
  float filterSum = 0.f;
  constint halfWidth = blurFilterWidth/2;
  for (int r = -halfWidth; r <= halfWidth; ++r) 
  {
    for (int c = -halfWidth; c <= halfWidth; ++c) 
    {
      float filterValue = expf( -(float)(c * c + r * r) / (2.f * blurFilterSigma * blurFilterSigma));
      (*filter)[(r + halfWidth) * blurFilterWidth + c + halfWidth] = filterValue;
      filterSum += filterValue;
    }
  }
  float normalizationFactor = 1.f / filterSum;
  for (int r = -halfWidth; r <= halfWidth; ++r) 
  {
    for (int c = -halfWidth; c <= halfWidth; ++c) 
    {
      (*filter)[(r + halfWidth) * blurFilterWidth + c + halfWidth] *= normalizationFactor;
    }
  }  
}
voidfreeFilter(float *filter){
  delete[] filter;
}
intmain( int argc, char** argv ){
  usingnamespace cv;
  usingnamespacestd;  
  usingnamespacestd::chrono;
  if( argc != 2)
  {
    cout <<" Usage: blur_image imagefile" << endl;
    return-1;
  }
  Mat image, blurredImage, referenceBlurredImage;
  uchar4 *imageArray, *blurredImageArray;
  prepareImagePointers(argv[1], image, &imageArray, blurredImage, &blurredImageArray, CV_8UC4);
  int numRows = image.rows, numCols = image.cols;
  float *filter, filterSigma;
  int filterWidth;
  prepareFilter(&filter, &filterWidth, &filterSigma);
  cv::Size filterSize(filterWidth, filterWidth);
  auto start = system_clock::now();
  cv::GaussianBlur(image, referenceBlurredImage, filterSize, filterSigma, filterSigma, BORDER_REPLICATE);
  auto duration = duration_cast<milliseconds>(system_clock::now() - start);
  cout<<"OpenCV time (ms):" << duration.count() << endl; 
  start = system_clock::now();
  BlurImageOpenMP(imageArray, blurredImageArray, numRows, numCols, filter, filterWidth);
  duration = duration_cast<milliseconds>(system_clock::now() - start);
  cout<<"OpenMP time (ms):" << duration.count() << endl; 
  cout<<"OpenMP similarity:" << getEuclidianSimilarity(referenceBlurredImage, blurredImage) << endl;
  for (int i=0; i<4; ++i)
  {
    memset(blurredImageArray, 0, sizeof(uchar4)*numRows*numCols); 
    start = system_clock::now();
    BlurImageCUDA(imageArray, blurredImageArray, numRows, numCols, filter, filterWidth);
    duration = duration_cast<milliseconds>(system_clock::now() - start);
    cout<<"CUDA time full (ms):" << duration.count() << endl;
    cout<<"CUDA similarity:" << getEuclidianSimilarity(referenceBlurredImage, blurredImage) << endl;
  }
  freeFilter(filter);
  return0;
}

1. We read the image file, prepare pointers to the original image and the resulting blurred image. The prepareImagePointers function remains the same, if necessary, you can see its source code on bitbucket.
2. We are preparing a Gaussian filter - that is, a set of our scales. We also remember the used filter parameters, so that later they are transferred to OpenCV and get a sample of a blurred image to check the correct operation of our algorithms.
3. We call the Gaussian blur function from OpenCV, save the resulting sample, measure the time spent.
4. We call the Gaussian blur function written using OpenMP, measure the time spent, compare the result with the sample. The image similarity calculation function getEuclidianSimilarity is as follows:
   
   getEuclidianSimilarity

   doublegetEuclidianSimilarity(const cv::Mat& a, const cv::Mat& b){
     double errorL2 = cv::norm(a, b, cv::NORM_L2);
     double similarity = errorL2 / (double) (a.rows * a.cols);
     return similarity;
   }
   

   
   
   
   In fact, it finds the average sum of squares of the differences in the values ​​of all channels of all pixels of two images.
5. We call the CUDA version of the Gaussian blur 4 times, each time measuring the time spent and comparing the result with the sample. Why call 4 times? The fact is that during the very first call, a certain time will be spent on initialization - therefore, it is better to start several times and measure the time spent on subsequent calls.

#include<stdio.h>#include<omp.h>#include<vector_types.h>#include<vector_functions.h>voidBlurImageOpenMP(const uchar4 * const imageArray, 
					 uchar4 * const blurredImageArray, 
					 constlong numRows, 
					 constlong numCols, 
					 constfloat * const filter, 
					 constsize_t filterWidth){
	usingnamespacestd;
	constlong halfWidth = filterWidth/2;
    #pragma omp parallel for collapse(2)for (long row = 0; row < numRows; ++row)
    {
        for (long col = 0; col < numCols; ++col)
        {
    	  float resR=0.0f, resG=0.0f, resB=0.0f;
		  for (long filterRow = -halfWidth; filterRow <= halfWidth; ++filterRow) 
		  {
		    for (long filterCol = -halfWidth; filterCol <= halfWidth; ++filterCol) 
		    {
		      //Find the global image position for this filter position//clamp to boundary of the imageconstlong imageRow = min(max(row + filterRow, static_cast<long>(0)), numRows - 1);
		      constlong imageCol = min(max(col + filterCol, static_cast<long>(0)), numCols - 1);
		      const uchar4 imagePixel = imageArray[imageRow*numCols+imageCol];
		      constfloat filterValue = filter[(filterRow+halfWidth)*filterWidth+filterCol+halfWidth];
		      resR += imagePixel.x*filterValue;
		      resG += imagePixel.y*filterValue;
		      resB += imagePixel.z*filterValue;
		    }
		  }
		  blurredImageArray[row*numCols+col] = make_uchar4(resR, resG, resB, 255);
        }
    }
}

#include<iostream>#include<cuda.h>#include"CUDA_wrappers.hpp"#include"common/CUDA_common.hpp"
__global__
voidgaussian_blur(const uchar4* const d_image,
                   uchar4* const d_blurredImage,
                   constint numRows,
                   constint numCols,
                   constfloat * const d_filter, 
                   constint filterWidth){
  constint row = blockIdx.y*blockDim.y+threadIdx.y;
  constint col = blockIdx.x*blockDim.x+threadIdx.x;
  if (col >= numCols || row >= numRows)
    return;
  constint halfWidth = filterWidth/2;
  extern __shared__ float shared_filter[]; 
  if (threadIdx.y < filterWidth && threadIdx.x < filterWidth)
  {
    constint filterOff = threadIdx.y*filterWidth+threadIdx.x;
    shared_filter[filterOff] = d_filter[filterOff];    
  }
  __syncthreads();
  float resR=0.0f, resG=0.0f, resB=0.0f;
  for (int filterRow = -halfWidth; filterRow <= halfWidth; ++filterRow) 
  {
    for (int filterCol = -halfWidth; filterCol <= halfWidth; ++filterCol) 
    {
      //Find the global image position for this filter position//clamp to boundary of the imageconstint imageRow = min(max(row + filterRow, 0), numRows - 1);
      constint imageCol = min(max(col + filterCol, 0), numCols - 1);
      const uchar4 imagePixel = d_image[imageRow*numCols+imageCol];
      constfloat filterValue = shared_filter[(filterRow+halfWidth)*filterWidth+filterCol+halfWidth];
      resR += imagePixel.x * filterValue;
      resG += imagePixel.y * filterValue;
      resB += imagePixel.z * filterValue;        
    }
  }
  d_blurredImage[row*numCols+col] = make_uchar4(resR, resG, resB, 255);     
}
voidBlurImageCUDA(const uchar4 * const h_image, 
                   uchar4 * const h_blurredImage, 
                   constsize_t numRows, 
                   constsize_t numCols, 
                   constfloat * const h_filter, 
                   constsize_t filterWidth){
  uchar4 *d_image, *d_blurredImage;
  cudaSetDevice(0);
  checkCudaErrors(cudaGetLastError());
  constsize_t numPixels = numRows * numCols;
  constsize_t imageSize = sizeof(uchar4) * numPixels;
  //allocate memory on the device for both input and output
  checkCudaErrors(cudaMalloc(&d_image, imageSize));
  checkCudaErrors(cudaMalloc(&d_blurredImage, imageSize));
  //copy input array to the GPU
  checkCudaErrors(cudaMemcpy(d_image, h_image, imageSize, cudaMemcpyHostToDevice));
  float *d_filter;
  constsize_t filterSize = sizeof(float) * filterWidth * filterWidth;
  checkCudaErrors(cudaMalloc(&d_filter, filterSize));
  checkCudaErrors(cudaMemcpy(d_filter, h_filter, filterSize, cudaMemcpyHostToDevice));
  dim3 blockSize;
  dim3 gridSize;
  int threadNum;
  cudaEvent_t start, stop;
  cudaEventCreate(&start);
  cudaEventCreate(&stop);
  threadNum = 32;
  blockSize = dim3(threadNum, threadNum, 1);
  gridSize = dim3(numCols/threadNum+1, numRows/threadNum+1, 1);
  cudaEventRecord(start);
  gaussian_blur<<<gridSize, blockSize, filterSize>>>(d_image, 
                                                     d_blurredImage, 
                                                     numRows, 
                                                     numCols, 
                                                     d_filter, 
                                                     filterWidth);
  cudaEventRecord(stop);
  cudaEventSynchronize(stop);
  cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
  float milliseconds = 0;
  cudaEventElapsedTime(&milliseconds, start, stop);
  std::cout << "CUDA time kernel (ms): " << milliseconds << std::endl;
  checkCudaErrors(cudaMemcpy(h_blurredImage, d_blurredImage, sizeof(uchar4) * numPixels, cudaMemcpyDeviceToHost));
  checkCudaErrors(cudaFree(d_filter));
  checkCudaErrors(cudaFree(d_image));
  checkCudaErrors(cudaFree(d_blurredImage));
}

1. Allocate memory on the device for the source image, output image, filter. Copy the corresponding data from the host memory to the allocated device memory.
2. We call the kernel. Pay attention to the new 3rd parameter when invoking the kernel: <<< gridSize, blockSize, filterSize >>> - it sets the size of the shared memory that each block needs. In this task, it was possible to implement the use of the block’s shared memory in two ways: either move the filter data to the block’s shared memory, or move there a “piece” of the image that is needed only for this block, since this data is needed by several block flows at once. However, the second option is somewhat more complicated - you need to consider that the piece of the input image that is needed for each block is somewhat larger than the block itself - after all, we perform the gather operation, which means that each stream calculates the values ​​of one pixel of the output image using several neighboring pixels of the original image:
   
   Therefore, I settled on the first option, which means that each block needs exactly sizeof (float) * filterWidth * filterWidth memory to store all filter values. Moving the filter weights from the device’s memory to the block’s shared memory occurs as follows:
   
   Hidden text

   extern __shared__ float shared_filter[]; 
     if (threadIdx.y < filterWidth && threadIdx.x < filterWidth)
     {
       constint filterOff = threadIdx.y*filterWidth+threadIdx.x;
       shared_filter[filterOff] = d_filter[filterOff];    
     }
     __syncthreads();
   

   
   
   
   Here __shared__ in the declaration of the array of filter weights says that this data should be placed in the shared memory of the block; extern means that the size of the allocated memory will be set at the kernel call point; __syncthreads is a barrier that guarantees that all filter weights will be transferred to the shared memory of a block before any of the threads of this block continues to execute. Further, all readings of the filter weights are carried out already from the general memory of the block.
3. Copy the output image from the device’s memory to the host’s memory, free the allocated memory.

OpenCV time (ms):2
OpenMP time (ms):11
OpenMP similarity:0.00287131
CUDA time kernel (ms): 2.93245
CUDA time full (ms):32
CUDA similarity:0.00287131
CUDA time kernel (ms): 2.93402
CUDA time full (ms):4
CUDA similarity:0.00287131
CUDA time kernel (ms): 2.93267
CUDA time full (ms):4
CUDA similarity:0.00287131
CUDA time kernel (ms): 2.93312
CUDA time full (ms):4
CUDA similarity:0.00287131