OneApi实验

OneAPI介绍

OneAPI是一种通用的并行编程库，可以适配Intel的多种异构硬件。这个库的主要用途是在高性能计算（HPC）场景以及深度学习的卷积优化、向量传播等场景中，降低适配不同硬件的编译适配成本。OneAPI还提供了一个云平台DevCloud，用户可以免费注册后使用，即使你的设备不是Intel的，也可以体验到Intel OneAPI配套硬件的性能，而且无需自行安装toolkit。我在本次实验中使用的云平台链接是jupyter.oneapi.devcloud.intel.com，本次实验主要基于Jupyter Lab进行，使用Jupyter Notebook构建了可视化界面，即使你不懂如何使用命令行进行编译，也可以直接复制官方教程中的脚本进行编译。

实验过程

本次实验的目标是计算矩阵乘法，需要编写一个基于oneAPI的C++/SYCL程序来执行矩阵乘法操作。这个问题的难点在于需要处理大尺寸矩阵的乘法操作以及不同线程之间的数据依赖关系。在实现矩阵乘法时，我选择使用块矩阵乘法以及共享内存来提高计算效率。
虽然有一些优化矩阵乘法的方法，例如Strassen矩阵乘法，但是这些方法往往有依赖关系，子问题数量也呈指数级增加，而我们的机器的线程数（或者说核数）是有限的，所以不能无限制的增加，将子问题派发给不同线程。因此，我选择使用复杂度为O(n^3)的GEMM矩阵乘法算法，通过并行计算来缩短计算的时间，也就是达成一个加速比。
在实现过程中，我选择对矩阵进行分块计算。这里只考虑了矩阵维度能对齐的情况，如果不能对齐tile（例如大小不是512，而是511，513，那么需要进行zero-padding）。实际上的分块（tile）大小和机器的缓存大小有关，因为需要考虑矩阵乘法的访问序列是一行与一列共同访问的，所以一旦tile太大，访问一列时就会频繁发生cache miss，访问A矩阵和B矩阵的行和列时也会频繁触发switch，影响性能。

代码

#include <chrono>
#include <iostream>
#include <sycl/sycl.hpp>

#define random_float() (rand() / double(RAND_MAX))

using namespace sycl;

const int tileX = 8;
const int tileY = 8;

double execute_parallel_matrix_multiplication(float *A, float *B, float *C, 
                  int M, int N, int K, 
                  int BLOCK, sycl::queue &q) {

  auto grid_rows = M / tileY;
  auto grid_cols = N / tileX;
  auto local_ndrange  = range<2>(BLOCK, BLOCK);
  auto global_ndrange = range<2>(grid_rows, grid_cols);

  double duration = 0.0f;

  auto e = q.submit([&](sycl::handler &h) {
      h.parallel_for(
          sycl::nd_range<2>(global_ndrange, local_ndrange), [=](sycl::nd_item<2> index) {

              int row = tileY * index.get_global_id(0);
              int col = tileX * index.get_global_id(1);

              float sum[tileY][tileX] = {0.0f};
              float subA[tileY] = {0.0f};
              float subB[tileX] = {0.0f};

              for (int k = 0; k < N; k++) {

                for(int m = 0; m < tileY; m++) {
                    subA[m] = A[(row + m) * N + k];
                } 

                for(int p = 0; p < tileX; p++) {
                    subB[p] = B[k * N + p + col];
                } 

                for (int m = 0; m < tileY; m++) {
                  for (int p = 0; p < tileX; p++) {
                    sum[m][p] += subA[m] * subB[p];
                  }
                }
              }

              for (int m = 0; m < tileY; m++) {
                for (int p = 0; p < tileX; p++) {
                  C[(row + m) * N + col + p] = sum[m][p];
                }
              }

          });
    });
    e.wait();

    duration += (e.get_profiling_info<info::event_profiling::command_end>() -
    e.get_profiling_info<info::event_profiling::command_start>()) /1000.0f/1000.0f;

    return(duration);
}

double execute_normal_matrix_multiplication(float *cA, float *cB, float *cC, int M, int N, int K) {
    
    double duration = 0.0;
    std::chrono::high_resolution_clock::time_point s, e;
    s = std::chrono::high_resolution_clock::now();
    for(int i = 0; i < M; i++) {
        for(int j = 0; j < N; j++) {
            float sum = 0.0f;
            for(int k = 0; k < K; k++) {
                sum +=  cA[i * K + k] * cB[k * N  + j];
            }
            cC[i * N + j] = sum;
        }
    }
    e = std::chrono::high_resolution_clock::now();
    duration = std::chrono::duration<float, std::milli>(e - s).count();
    return(duration);
}

int verify_results(float *normal_res, float *parallel_res, int length){
    int err = 0;
    for(int i = 0; i < length; i++) {
       if( fabs(normal_res[i] - parallel_res[i]) > 1e-3) {
          err++;
          printf("\n%lf, %lf, %d %lf", normal_res[i], parallel_res[i], i, fabs(normal_res[i]-parallel_res[i]));
       } 
    }
    return(err);
}

int execute_gemm(const int M, 
         const int N, 
         const int K, 
         const int block_size,
         const int iterations, 
         sycl::queue &q) {

  std::cout << "Problem size: c(" << M << "," <<  N << ") ="
       << " a(" << M << "," << K << ") *" 
       << " b(" << K << "," << N << ")\n";

  auto A = malloc_shared<float>(M * K, q);
  auto B = malloc_shared<float>(K * N, q);
  auto C = malloc_shared<float>(M * N, q);
  auto C_host = malloc_host<float>(M * N, q);

  for(int i=0; i < M * K; i++) {
      A[i] = random_float();
  }

  for(int i=0; i < K * N; i++) {
      B[i] = random_float();
  }

  for(int i=0; i < M * N; i++) {
      C[i] = 0.0f;
      C_host[i] = 0.0f;
  }

  double flopsPerMatrixMul
      = 2.0 * static_cast<double>(M) * static_cast<double>(N) * static_cast<double>(K);

  double duration_parallel = 0.0f;
  double duration_normal = 0.0f;

  int warmup = 10;
  for (int run = 0; run < iterations + warmup; run++) {
    float duration = execute_parallel_matrix_multiplication(A, B, C, M, N, K, block_size, q);
    if(run >= warmup) duration_parallel += duration;
  }
  duration_parallel = duration_parallel / iterations;

  warmup = 2;
  for(int run = 0; run < iterations/2 + warmup; run++) {
      float duration = execute_normal_matrix_multiplication(A, B, C_host, M, N, K);
      if(run >= warmup) duration_normal += duration;
  }
  duration_normal = duration_normal / iterations/2;

  int errCode = 0;
  if(errCode > 0) printf("\nThere are %d errors\n", errCode);

  printf("\nGEMM size M = %d, N = %d, K = %d", M, N, K);
  printf("\nWork-Group size = %d * %d, tile_X = %d, tile_Y = %d", block_size, block_size, tileX, tileY);
  printf("\nPerformance Flops = %lf, \n" 
          "Parallel Computation Time = %lf (ms); \n"
          "Normal Computaiton Time = %lf (ms); \n"
          "Speedup = %lf\n", 
          flopsPerMatrixMul, duration_parallel, duration_normal, duration_normal/duration_parallel);

  free(A, q);
  free(B, q);
  free(C, q);
  free(C_host, q);

  return(errCode);
}

int main() {
    
  auto propList = sycl::property_list {sycl::property::queue::enable_profiling()};
  queue my_queue(default_selector_v , propList);

  int errCode = execute_gemm(512, 512, 512, 
                     4,             
                     10,               
                     my_queue);

  return(errCode);
}

总结

通过采用Intel® oneAPI Data Parallel C++ (DPC++) 编译器来优化矩阵乘法运算，可以实现显著的性能提升。这种加速对于那些需要处理巨量数据集的应用来说至关重要，因为它不仅大幅度缩短了计算时间，同时还提高了能效，使得数据分析、机器学习算法训练以及复杂科学模拟等计算密集型任务变得更加高效和可行。此外，利用Intel® oneAPI DPC++编译器的跨平台特性，开发者能够针对多种架构（包括CPU、GPU和FPGA）编写统一的代码，从而进一步提升软件的可移植性和灵活性。