0
点赞
收藏
分享

微信扫一扫

cuda编程入门（3）—— 实现矩阵乘法(一)

蚁族的乐土 2022-06-28 阅读 73
矩阵算法线性代数矩阵乘法#includeOpenStack云计算


文章目录

  • ​​一、CPU下的一般矩阵乘法​​
  • ​​二、CPU下循环交换矩阵乘法​​
  • ​​三、CPU下转置矩阵乘法​​
  • ​​附录（连续访问为什么比随机访问快）​​

本文章为 ​​《 GPU编程与优化 大众高性能计算》​​ 的读书笔记，例子也都取自书中教程。

一、CPU下的一般矩阵乘法

cuda编程入门（3）—— 实现矩阵乘法(一)_算法

// 矩阵乘法一
// m*l l*n
void matrix_mul(float* x, float * y, float* z, int m, int n, int l)
{
  for(int i = 0; i < m; i++) {
    for(int j = 0; j < n; j++) {
      for(int k = 0; k < l; k++) {
        z[i*n+ j] += x[i*l + k] * y[k*n + j];
      }
    }
  }
}

#include <iostream>
#include <assert.h>
#include <chrono>

int main()
{
  int M = 2048;
  int L = 1024;
  int N = 512;

  // 申请host内存
  float *x = NULL;
  float *y = NULL;
  float *z = NULL;
  x = (float*)malloc(M*L*sizeof(float));
  y = (float*)malloc(L*N*sizeof(float));
  z = (float*)malloc(M*N*sizeof(float));

  if(x == NULL || y == NULL || z == NULL)
    return 0;

  // 初始化数据
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < L; ++j) {
      x[i*L + j] = 1.1;
    }
  }
  for (int i = 0; i < L; ++i) {
    for (int j = 0; j < N; ++j) {
      y[i*N + j] = 1.1;
    }
  }
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < N; ++j) {
      z[i*N + j] = 0;
    }
  }

  // 执行矩阵乘法
  using namespace std::chrono;
  steady_clock::time_point t1 = steady_clock::now();
  matrix_mul(x, y, z, M, N, L);
  steady_clock::time_point t2 = steady_clock::now();

  duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
  std::cout << "It took me " << time_span.count() << " seconds.";
  std::cout << std::endl;

  // 输出前10个数值
  for(int i = 0; i < 10; i++) {
    std::cout << z[i] << " ";
  }
  std::cout << std::endl;

  std::cout << "Done!" << std::endl;

  // 释放host内存
  free(x);
  free(y);
  free(z);

  return 0;
}

matrix_mul函数中 对矩阵 y 的数据访问是不连续的，会影响速度。

二、CPU下循环交换矩阵乘法

cuda编程入门（3）—— 实现矩阵乘法(一)_矩阵乘法_02

// 矩阵乘法一
// m*l l*n
void matrix_mul(float* x, float * y, float* z, int m, int n, int l)
{
  for(int i = 0; i < m; i++) {
    for(int k = 0; k < l; k++) {
      for(int j = 0; j < n; j++) {
        z[i*n + j] += x[i*l + k] * y[k*n + j];
      }
    }
  }
}

#include <iostream>
#include <assert.h>
#include <chrono>

int main()
{
  int M = 2048;
  int L = 1024;
  int N = 512;

  // 申请host内存
  float *x = NULL;
  float *y = NULL;
  float *z = NULL;
  x = (float*)malloc(M*L*sizeof(float));
  y = (float*)malloc(L*N*sizeof(float));
  z = (float*)malloc(M*N*sizeof(float));

  if(x == NULL || y == NULL || z == NULL)
    return 0;

  // 初始化数据
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < L; ++j) {
      x[i*L + j] = 1.1;
    }
  }
  for (int i = 0; i < L; ++i) {
    for (int j = 0; j < N; ++j) {
      y[i*N + j] = 1.1;
    }
  }
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < N; ++j) {
      z[i*N + j] = 0;
    }
  }

  // 执行矩阵乘法
  using namespace std::chrono;
  steady_clock::time_point t1 = steady_clock::now();
  matrix_mul(x, y, z, M, N, L);
  steady_clock::time_point t2 = steady_clock::now();

  duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
  std::cout << "It took me " << time_span.count() << " seconds.";
  std::cout << std::endl;

  // 输出前10个数值
  for(int i = 0; i < 10; i++) {
    std::cout << z[i] << " ";
  }
  std::cout << std::endl;

  std::cout << "Done!" << std::endl;

  // 释放host内存
  free(x);
  free(y);
  free(z);

  return 0;
}

三、CPU下转置矩阵乘法

cuda编程入门（3）—— 实现矩阵乘法(一)_算法_03

// 矩阵乘法一
// m*l l*n
#include <stdlib.h>
void matrix_mul(float* x, float * y, float* z, int m, int n, int l)
{
  float* transposed = (float*)malloc(l*n*sizeof(float));
  // l*n --> n*l
  for(int i = 0; i < n; i++) {
    for(int j = 0; j < l; j++) {
      transposed[i*l + j] = y[j*n + i];
    }
  }

  //m*l n*l
  for(int i = 0; i < m; i++) {
    for(int j = 0; j < n; j++) {
      for(int k = 0; k < l; k++) {
        z[i*n + j] += x[i*l + k] * y[j*l + k];
      }
    }
  }

  free(transposed);
}

#include <iostream>
#include <assert.h>
#include <chrono>

int main()
{
  int M = 2048;
  int L = 1024;
  int N = 512;

  // 申请host内存
  float *x = NULL;
  float *y = NULL;
  float *z = NULL;
  x = (float*)malloc(M*L*sizeof(float));
  y = (float*)malloc(L*N*sizeof(float));
  z = (float*)malloc(M*N*sizeof(float));

  if(x == NULL || y == NULL || z == NULL)
    return 0;

  // 初始化数据
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < L; ++j) {
      x[i*L + j] = 1.1;
    }
  }
  for (int i = 0; i < L; ++i) {
    for (int j = 0; j < N; ++j) {
      y[i*N + j] = 1.1;
    }
  }
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < N; ++j) {
      z[i*N + j] = 0;
    }
  }

  // 执行矩阵乘法
  using namespace std::chrono;
  steady_clock::time_point t1 = steady_clock::now();
  matrix_mul(x, y, z, M, N, L);
  steady_clock::time_point t2 = steady_clock::now();

  duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
  std::cout << "It took me " << time_span.count() << " seconds.";
  std::cout << std::endl;

  // 输出前10个数值
  for(int i = 0; i < 10; i++) {
    std::cout << z[i] << " ";
  }
  std::cout << std::endl;

  std::cout << "Done!" << std::endl;

  // 释放host内存
  free(x);
  free(y);
  free(z);

  return 0;
}

附录（连续访问为什么比随机访问快）

​​内存（DRAM）的连续读写速度和随机读写速度是一样的吗？​​

作者：独孤星夜
链接：https://www.zhihu.com/question/325168076/answer/689074506
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

内存每一个bank只有一个buffer，读写都要通过这个buffer来进行，这个buffer只能装载内存的一行，也叫page。存在buffer中的行的状态称为page opened，而同一bank其他行的状态称为page closed。如果到来的命令正好是访问当前open状态的行，那么这叫page hit，只需要一个cas命令(列选择，既读或者写）即可完成操作。如果所有行处于close状态，这叫page miss，需要先发一个ras命令，也叫active命令，行选择。然后再发cas命令完成操作。如果当前访问的行和open状态的行不是同一行，那么这叫page conflict。需要先发一个pre-charge命令，关闭当前open的行，然后再ras命令，最后cas命令完成操作。所以memory controller主要使用两种page策略，open page 和close page。open page，顾明思议，读写后不会立刻进行pre-charge，这种策略对连续地址读写方式更友好，但一旦出现page conflict，受到的影响就很明显。而close page则是每次读写后立刻pre-charge，这样每次到来的命令都是page miss，不会出现page conflict，对随机读写更加友好。一般来说，连续命令以连续地址操作较多，所以一般性能好优化好的memory controller采用open page方式的占多数，需要对命令进行重排序，再调度，保持同一地址的读写顺序，pre-charge预测等等优化，还要采取对应的地址映射方式（这块挺复杂，简单的说，可以把collum地址映射到地址的低位，采用CS interleave，把bank group和bank地址映射到低位，bank group比bank更低为好），目的是增加page hit，减少page conflict。close page方式是结构相对简单的memory controller主要采用。所以连续地址读写和随机地址读写的速度，根据memory controller采用的page策略有关，当然大部分是连续地址读写速度更快带宽更高。另外，从读命令到写命令所用切换时间比从写到读少，这是因为数据操作其实主要在buffer中进行，如果前一个是读，读命令不会影响数据，buffer中数据始终和open的行中数据保持一致，所以pre-charge的时候不需要把buffer中数据写回行中。而写命令会改变数据，需要一个写回的过程。所以还尽量要减少写到读的切换，在保证正确读写顺序的基础上。


举报

相关推荐

cuda编程入门（4）—— 实现矩阵乘法(二)

矩阵算法c++矩阵乘法i++虚拟化云计算

矩阵乘法实现

算法

[8] CUDA之向量点乘和矩阵乘法

人工智能搜索引擎

cuda编程入门（1）

深度学习人工智能多线程核函数数据Java编程语言

hdu4920 矩阵乘法%3

#include#pragma时间复杂度编程语言

矩阵乘法

矩阵乘法i++快速幂斐波那契数列编程语言

cuda从入门到精通（六）共享内存和循环分块实现CUDA矩阵乘

矩阵算法性能优化

python矩阵乘法

python

矩阵乘法(java)

i++Java编程语言

pytorch矩阵乘法

linux
0 条评论