基于DPC++编写矩阵乘法计算代码

最新推荐文章于 2024-09-16 17:46:36 发布
最新推荐文章于 2024-09-16 17:46:36 发布
阅读量78 收藏 1
点赞数
文章标签： c++
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/m0_61726571/article/details/131272465
版权
该代码示例展示了使用CUDA进行矩阵乘法时，通过调整分块尺寸(tile_X,tile_Y)来优化GPU计算性能的过程。通过对比不同尺寸的分块，发现在tile_X=8,tile_Y=8时，计算速度达到最优。这表明正方形分块在某些情况下能提供更好的并行计算效率。程序通过比较GPU和CPU的计算时间，并验证了结果的一致性，以评估优化效果。
摘要由CSDN通过智能技术生成

Basic_parafor示例

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

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

using namespace std;
using namespace sycl;

// return execution time
double gpu_kernel(float *A, float *B, float *C, int M, int N, int K, int block_size, sycl::queue &q) {

  // define the workgroup size and mapping
  auto grid_rows = (M + block_size - 1) / block_size * block_size;
  auto grid_cols = (N + block_size - 1) / block_size * block_size;
  auto local_ndrange  = range<2>(block_size, block_size);
  auto global_ndrange = range<2>(grid_rows, grid_cols);

  double duration = 0.0f;
  auto e = q.submit([&](sycl::handler &h) {
      h.parallel_for<class k_name_t>(
          sycl::nd_range<2>(global_ndrange, local_ndrange), [=](sycl::nd_item<2> index) {

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

              float sum = 0.0f;

              for (int i = 0; i < K; i++) {
                sum += A[row * K + i] * B[i * N  + col];
              }
              C[row * N + col] = sum;
          });
    });
    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);
}

// return execution time
double cpu_kernel(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;

    // Single Thread Computation in CPU 
    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(float *cpu_res, float *gpu_res, int length){
    int err = 0;
    for(int i = 0; i < length; i++) {
       if( fabs(cpu_res[i] - gpu_res[i]) > 1e-3) {
          err++;
          printf("\n%lf, %lf", cpu_res[i], gpu_res[i]);
       } 
    }
    return(err);
}

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

  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);

  // init the A/B/C
  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_gpu = 0.0f;
  double duration_cpu = 0.0f;

  // GPU compuation and timer 
  int warmup = 10;
  for (int run = 0; run < iterations + warmup; run++) {
    float duration = gpu_kernel(A, B, C, M, N, K, block_size, q);
    if(run >= warmup) duration_gpu += duration;
  }
  duration_gpu = duration_gpu / iterations;

  // CPU compuation and timer 
  warmup = 2;
  for(int run = 0; run < iterations/2 + warmup; run++) {
      float duration = cpu_kernel(A, B, C_host, M, N, K);
      if(run >= warmup) duration_cpu += duration;
  }
  duration_cpu = duration_cpu / iterations/2;

  // Compare the resutls of CPU and GPU 
  int errCode = 0;
  errCode = verify(C_host, C, M*N);
  if(errCode > 0) printf("\nThere are %d errors\n", errCode);

  printf("\nPerformance Flops = %lf, \n" 
          "GPU Computation Time = %lf (ms); \n"
          "CPU Computaiton Time = %lf (ms); \n", 
          flopsPerMatrixMul, duration_gpu, duration_cpu);

  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_gpu_queue( sycl::gpu_selector{} , propList);

  int errCode = gemm(2000, 2000, 2000, 4, 10, my_gpu_queue);

  return(errCode);
}

修改程序输入数据的大小，设定M=N=K=2000

测试不同tile_X和tile_Y大小下，矩阵计算的性能


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

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

using namespace std;
using namespace sycl;

#define tileY 8
#define tileX 8

// return execution time
double gpu_kernel(float *A, float *B, float *C, 
                  int M, int N, int K, 
                  int BLOCK, sycl::queue &q) {

  // define the workgroup size and mapping
  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<class k_name_t>(
          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};

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

                // read data to register
                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];
                  }
                }

              } //end of K

              // write results back
              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);
}

// return execution time
double cpu_kernel(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;

    // Single Thread Computation in CPU 
    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(float *cpu_res, float *gpu_res, int length){
    int err = 0;
    for(int i = 0; i < length; i++) {
       if( fabs(cpu_res[i] - gpu_res[i]) > 1e-3) {
          err++;
          printf("\n%lf, %lf", cpu_res[i], gpu_res[i]);
       } 
    }
    return(err);
}

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

  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);

  // init the A/B/C
  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_gpu = 0.0f;
  double duration_cpu = 0.0f;

  // GPU compuation and timer 
  int warmup = 10;
  for (int run = 0; run < iterations + warmup; run++) {
    float duration = gpu_kernel(A, B, C, M, N, K, block_size, q);
    if(run >= warmup) duration_gpu += duration;
  }
  duration_gpu = duration_gpu / iterations;

  // CPU compuation and timer 
  warmup = 2;
  for(int run = 0; run < iterations/2 + warmup; run++) {
      float duration = cpu_kernel(A, B, C_host, M, N, K);
      if(run >= warmup) duration_cpu += duration;
  }
  duration_cpu = duration_cpu / iterations/2;

  // Compare the resutls of CPU and GPU 
  int errCode = 0;
  errCode = verify(C_host, C, M*N);
  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" 
          "GPU Computation Time = %lf (ms); \n"
          "CPU Computaiton Time = %lf (ms); \n", 
          flopsPerMatrixMul, duration_gpu, duration_cpu);

  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_gpu_queue( sycl::gpu_selector_v , propList);

  int errCode = gemm(512, 512, 512, /* GEMM size, M, N, K */
                     4,             /* workgroup size */ 
                     10,            /* repeat time */   
                     my_gpu_queue);

  return(errCode);
}
  1. tile_X=2, tile_Y=2

  1. tile_X=8 tile_Y=8

  1. tile_X=16, tile_Y=16

  1. tile_X=16, tile_Y=8

  1. tile_X=4, tile_Y=8

由上述结果，当tile_X=8 tile_Y=8时，计算速度最优，分块尺寸为正方形的速度快与分块尺寸为非方形的速度

 

 

m0_61726571
关注
矩阵链乘法C++实现
OOM
05-08 2163
#include #define MAX_SIZE 100 using namespace std; int num; int size[MAX_SIZE],s[MAX_SIZE][MAX_SIZE]; int ji[MAX_SIZE][MAX_SIZE]; void opt_matrix_muln(int size[]) { int n=num; int i,j,k,w,q
基于DPCPP编写矩阵乘法计算代码
wuwangwoc的博客
06-18 119
测试不同tile_X和tile_Y大小下（ tile_X和tile_Y可以不一致），矩阵计算的性能。2、对比（2，2），（2，4），（2，8）三组实验，y数值变大对于GPU的时间有影响但并不大。3、对比（2，2），（4，2），（8，2）三组实验，x数值变大对于GPU的时间影响很大。可通过修改MNK值（此例设置M=N=K=2000）决定矩阵参数，实现不同结果进行比较。编写矩阵乘法代码实现，并编译执行；对代码进行执行时间分析，比较不同实现的效率差异；1、x,y的值的变化对CPU时间影响不大。
c++矩阵连乘的动态规划算法并输出_算法面试必修课，动态规划基础题型归纳（三）
weixin_39960019的博客
11-28 401
动态规划（Dynamic Programming，简称DP），是大家都觉得比较难以掌握的算法。为了应付面试，我们经常会背诵一下DP问题的源码，其实，只要理解了思想，掌握基本的模型，然后再来点写代码的套路，动态规划并没有那么难。先来回顾一下题型目录。动态规划基本题型总结提纲目录1.硬币找零，路径规划2.字符串相似度/编辑距离（edit distance）3.最长公共子序列(Longest Commo...
【计算机算法设计与分析】矩阵连乘问题（C++_DP）
ccql's Blog
11-22 967
问题描述 给定n个矩阵｛A1,A2,…,An｝，其中Ai与Ai+1是可乘的，i=1，2…，n-1。如何确定计算矩阵连乘积的计算次序，使得依此次序计算矩阵连乘积需要的数乘次数最少。例如，给定三个连乘矩阵{A1，A2，A3}的维数分别是10100，1005和550，采用（A1A2）A3，乘法次数为101005+10550=7500次，而采用A1（A2A3），乘法次数为100550+10100*50=75000次乘法，显然，最好的次序是（A1A2)A3，乘法次数为7500次。 求解过程 Code #includ
矩阵链乘法 | Matrix Chain Multiplication | 动态规划 | C/C++实现
qq_40998706的博客
01-28 1654
问题描述 对于给定的n个矩阵形成的矩阵链M1,M2,M3,...,MnM_1,M_2,M_3,...,M_nM1​,M2​,M3​,...,Mn​，求计算乘积M1M2M3,...MnM_1M_2M_3,...M_nM1​M2​M3​,...Mn​时最少次标量相乘的运算顺序，这类问题称为矩阵链乘法问题(Matrix Chain Multiplication Problem)。 请编写一个程序，当给定...
矩阵链乘（C++实现）
clearlove_black的博客
05-27 1737
矩阵连乘 递推式 m[i][j] = max(m[i][j], m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j]) // 有点像分治法的思想，因为先分开最小的两个，之后不断合并成最大的一部分，也就是一个整体 #include &lt;iostream&gt; using namespace std; int m[100][100]; in...
通过DPC ++从CUDA到SYCL 2020进行光线追踪
04-11
本文将详细介绍如何将使用CUDA编写的光线追踪代码转换为DPC++，进而利用SYCL 2020的标准特性。首先，我们需要理解CUDA编程的基本概念，例如设备和主机的概念，以及kernel函数的定义和执行。CUDA代码通常包含host代码...
英特尔oneAPI——通过DPC++实现异构计算
weixin_44641025的博客
11-06 2187
OneAPI是英特尔提出的一种开放的跨架构的编程模型，使开发人员可以自由地跨多个架构使用单个代码库，从而帮助加速计算。 本文简要介绍了OneAPI和DPC++，并通过一个简单的实验介绍如何用DPC++实现异构并行计算程序的设计。
oneAPI DPC++并行计算的核函数示例
m0_68461595的博客
06-03 271
在这个示例中，我们使用 DPC++ 编写了一个名为 vectorAdd 的函数，该函数接受两个输入矢量 A 和 B，将它们相加并将结果存储在输出矢量 C 中。在并行任务中，我们使用 sycl::parallel_for 函数定义了一个内核函数，该函数使用索引 idx 并行地将 A 和 B 中的元素相加，并将结果存储在 C 中。通过 DPC++，开发人员可以使用统一的编程模型在不同的硬件加速器上实现高性能计算，包括在 CPU 和 GPU 上进行并行计算，以及在 FPGA 上进行定制化加速。
oneDPL:oneAPI DPC ++库（oneDPL）https
03-28
oneAPI DPC ++库（oneDPL）旨在与oneAPI DPC ++编译器一起使用，以向开发人员提供高生产率的API，从而可以最大程度地减少跨设备的DPC ++编程工作，以实现高性能并行应用程序。 oneDPL包含以下组件： 适用于DPC ++...
基于matlab+C/C++实现的K-means+FCM+谱聚类+DBSCAN+AP+DPC聚类算法比较+源码（期末大作业）
08-04
基于matlab+C/C++实现的K-means+FCM+谱聚类+DBSCAN+AP+DPC聚类算法比较+源码，适合毕业设计、课程设计、项目开发。项目源码已经过严格测试，可以放心参考并在此基础上延申使用，详情见md文档~ 基于matlab+C/C++实现...
DP—矩阵链乘法（代码）
pspdragon的博客
11-22 495
#include <iostream> int m[1003][1003]; //m[i][j]记录从第i个矩阵到第j个矩阵最少乘法执行次数 int s[1003][1003]; //s[i][j]记录从第i个矩阵到第j个矩阵最优质m[i][j]的分割点k int p[1003]; //第i个矩阵的行列分别是p[i-1]和p[i],1<=i<=n using namespace std; void
C++ STL中sort函数
weixin_43349440的博客
09-15 269
刷题中经常遇到，用于各种元素排序，可以自定义排序方式，vector、deque适用。详细讲解：https://zhuanlan.zhihu.com/p/36274119。，为避免QuickSort快排的递归调用带来过大的额外负荷，就改用。一旦分段后的数据量小于某个门槛。STL的sort算法，
CF C. Nikita and LCM
qq_74190237的博客
09-13 270
思维题.首先先对数组排序,考虑一种特殊情况,那就是如果数组最大值不能被其他数整除,这代表什么?代表其他数里面有某个质因子,最大值里面没有或者数量不管,所有最小公倍数都会大于最大值.那么其他情况的话,就可以确定所有的数都是最大数的因子,那么形成的最小公倍数也一定是最大数的因子,那么就可以枚举一下最大数的因子,就可以了.给长度为n的数组,求出最长的特殊子数组,如果子数组的最小公倍数不在原数组中,那么就是特殊的.
C++与QML的数据交互
qq_50679120的博客
09-11 526
这里没什么好说的，在C++中组合好json后，通过信号发送到qml中，qml中是可以解析json的。通过Q_PROPERTY属性，可以轻松地在QML中读取和修改数据，并将其更新到C++中，先创建一个类。在QML接收，直接通过data.can_id这样的形式去访问结构体中的元素。这两个也没什么好说的，在C++中直接发送信号就好了，在C++中声明并注册一个结构体。然后发送信号到QML中。
C++_map_set详解
最新发布
2302_77644537的博客
09-16 641
set是按照一定次序存储元素的容器在set中，元素的value也标识它(value就是key，类型为T)，并且每个value必须是唯一的。set中的元素不能在容器中修改(元素总是const)，但是可以从容器中插入或删除它们。在内部，set中的元素总是按照其内部比较对象(类型比较)所指示的特定严格的排序准则（默认是小于）进行排序。set在底层是用二叉搜索树(红黑树)实现的。map是关联容器，它按照特定的次序(按照key来比较)存储由键值key和值value组合而成的元素。
C语言从头学59——学习头文件math.h（二）
m0_72128260的博客
09-11 558
继续学习头文件的数学函数。
C++ 中 vector 的常用功能介绍
2202_75569688的博客
09-12 429
在 C++ 中，vector是一种常用的动态数组容器，提供了方便的自动扩展、内存管理以及各种便捷的操作方法。它是 C++ 标准模板库（STL）的一部分，适用于需要动态存储和管理大量元素的场景。在本文中，我们将简要介绍vector的常用功能，展示如何对其进行操作和排序。为了简洁起见，假设我们已经使用了。
DPC++ linux 环境变量
09-22
DPC在Linux环境中，需要设置DM_HOME环境变量和动态链接库文件的加载路径。可以通过以下步骤进行配置： 1. 切换至dmdba用户：`su - dmdba` 2. 执行命令使环境变量生效：`source .bash_profile` 3. 切换到root用户并...
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值