CUDA kernel to compare pairs of matrices

Question

My first time writing anything significant in CUDA.

This kernel takes two arrays representing square matrices and compares them pair-wise. It takes into consideration large input arrays, and calculates the appropriate number of threads on the GPU to spawn.

Currently takes ~1 second to compare 10 million 2x2 matrices (so 40 million floating point operations) on an RTX3060. This timing is inclusive of the time it takes to copy data to the GPU.

I feel this is a little slow for such a simple operation. Any suggestions to improve?

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>

// Kernel function to compare multiple pairs of square matrices element-wise
__global__ void compare_multiple_matrices(bool *results, const float *matricesA, const float *matricesB, int rank, int num_matrices, int matrices_offset)
{
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int matrix = blockIdx.z * blockDim.z + threadIdx.z + matrices_offset;

    results[matrix] = 1;

    if (row < rank && col < rank && matrix < num_matrices)
    {
        int index = matrix * rank * rank + row * rank + col;
        int is_equal = fabsf(matricesA[index] - matricesB[index]) <= 0.0001f;
        if (!is_equal)
        {
            results[matrix] = 0;
        }
    }
}

extern "C"
{
    int cuda_compare_multiple_matrices(bool *results, const float *matricesA, const float *matricesB, const int rank, const int num_matrices)
    {
        int res = 0;

        // Choose which GPU to run on, change this on a multi-GPU system.
        cudaError_t cudaStatus = cudaSetDevice(0);
        if (cudaStatus != cudaSuccess)
        {
            fprintf(stderr, "cudaSetDevice failed!  Do you have a CUDA-capable GPU installed?");

            return 1;
        }

        // Get the GPU device properties
        cudaDeviceProp deviceProp;
        int device;
        cudaGetDevice(&device);
        cudaGetDeviceProperties(&deviceProp, device);

        // Calculate the optimal blockDim.z based on the maximum thread count per block
        int blockDim_z = deviceProp.maxThreadsPerBlock / (16 * 16);

        dim3 blockDim(16, 16, blockDim_z);

        // Calculate the grid dimensions based on the GPU's maximum blocks per grid
        int gridDim_x = min((rank + blockDim.x - 1) / blockDim.x, deviceProp.maxGridSize[0]);
        int gridDim_y = min((rank + blockDim.y - 1) / blockDim.y, deviceProp.maxGridSize[1]);
        int gridDim_z = min((num_matrices + blockDim.z - 1) / blockDim.z, deviceProp.maxGridSize[2]);

        dim3 gridDim(gridDim_x, gridDim_y, gridDim_z);

        // Allocate GPU buffers for three vectors (two input, one output)    .
        float *dev_matricesA = 0;
        float *dev_matricesB = 0;
        bool *dev_results = 0;

        try
        {
            cudaStatus = cudaMalloc((void **)&dev_matricesA, num_matrices * rank * rank * sizeof(float));
            if (cudaStatus != cudaSuccess)
            {
                fprintf(stderr, "cudaMalloc failed!");
                throw;
            }

            cudaStatus = cudaMalloc((void **)&dev_matricesB, num_matrices * rank * rank * sizeof(float));
            if (cudaStatus != cudaSuccess)
            {
                fprintf(stderr, "cudaMalloc failed!");
                throw;
            }

            cudaStatus = cudaMalloc((void **)&dev_results, num_matrices * sizeof(bool));
            if (cudaStatus != cudaSuccess)
            {
                fprintf(stderr, "cudaMalloc failed!");
                throw;
            }

            // Copy input vectors from host memory to GPU buffers.
            cudaStatus = cudaMemcpy(dev_matricesA, matricesA, num_matrices * rank * rank * sizeof(float), cudaMemcpyHostToDevice);
            if (cudaStatus != cudaSuccess)
            {
                fprintf(stderr, "cudaMemcpy failed!");
                throw;
            }

            cudaStatus = cudaMemcpy(dev_matricesB, matricesB, num_matrices * rank * rank * sizeof(float), cudaMemcpyHostToDevice);
            if (cudaStatus != cudaSuccess)
            {
                fprintf(stderr, "cudaMemcpy failed!");
                throw;
            }

            int matrices_processed = 0;

            while (matrices_processed < num_matrices)
            {
                int matrices_remaining = num_matrices - matrices_processed;
                int matrices_to_process = min(matrices_remaining, blockDim.z * gridDim.z);

                compare_multiple_matrices<<<gridDim, blockDim>>>(dev_results, dev_matricesA, dev_matricesB, rank, num_matrices, matrices_processed);

                // Check for any errors launching the kernel
                cudaStatus = cudaGetLastError();
                if (cudaStatus != cudaSuccess)
                {
                    fprintf(stderr, "kernel launch failed: %s\n", cudaGetErrorString(cudaStatus));
                    throw;
                }

                // cudaDeviceSynchronize waits for the kernel to finish, and returns any errors encountered during the launch.
                cudaStatus = cudaDeviceSynchronize();
                if (cudaStatus != cudaSuccess)
                {
                    fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching kernel!\n", cudaStatus);
                    throw;
                }

                matrices_processed += matrices_to_process;
            }

            // Copy output vector from GPU buffer to host memory.
            cudaStatus = cudaMemcpy(results, dev_results, num_matrices * sizeof(bool), cudaMemcpyDeviceToHost);
            if (cudaStatus != cudaSuccess)
            {
                fprintf(stderr, "cudaMemcpy failed!");
                throw;
            }
        }
        catch (...)
        {
            res = 1;
        }

        cudaFree(dev_matricesA);
        cudaFree(dev_matricesB);
        cudaFree(dev_results);

        // cudaDeviceReset must be called before exiting in order for profiling and
        // tracing tools such as Nsight and Visual Profiler to show complete traces.
        cudaStatus = cudaDeviceReset();
        if (cudaStatus != cudaSuccess)
        {
            fprintf(stderr, "cudaDeviceReset failed!");

            res = 1;
        }

        return res;
    }
}
````

Answer 1

... this is a little slow for such a simple operation. Any suggestions to improve?

Avoid work

Minor idea: defer computations until needed with re-ordered code. Only assign once.

__global__ void compare_multiple_matrices(bool *results, const float *matricesA,
    const float *matricesB, int rank, int num_matrices, int matrices_offset) {
  int row = blockIdx.y * blockDim.y + threadIdx.y;
  if (row < rank) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    if (col < rank) {
      int matrix = blockIdx.z * blockDim.z + threadIdx.z + matrices_offset;
      if (matrix < num_matrices) {
        // int index = matrix * rank * rank + row * rank + col;
        int index = (matrix * rank + row) * rank + col;
        int is_equal = fabsf(matricesA[index] - matricesB[index]) <= 0.0001f;
        if (!is_equal) {
          results[matrix] = 0;
          return;
        }
      }
    }
  }
  results[matrix] = 1;
}

---

Non-speed issues

Beware overflow

Allocation uses: num_matrices * rank * rank * sizeof(float) which is ((int * int) * int) * size_t resulting in a size_t product. Avoid int overflow in the 1st 3 terms and use the widest type first.

// num_matrices * rank * rank * sizeof(float)
sizeof(float) * num_matrices * rank * rank

Likely not now a problem with 40 million, yet code may stumble with later/larger tasks.

Allocate to the referenced object, not type

I have found sizing the the reference object easier to code right, review and maintain and so reccomend:

// cudaStatus = cudaMalloc((void **)&dev_matricesA, 
//     num_matrices * rank * rank * sizeof(float));
cudaStatus = cudaMalloc((void **)&dev_matricesA, 
    sizeof dev_matricesA[0] * num_matrices * rank * rank);