code at: /lesson05/src_mat/src, data at: /lesson05/src_mat/data

Demo

// openmp_demo.cpp
#include <cstdio>
#include <omp.h>

int main(int argc, char **argv) {
	// 在这里使用这种指令开启 OpenMP
  #pragma omp parallel
  printf("Hello, OpenMP!\\n");
  return 0;
}

编译添加 OpenMP 相关库选项： g++ -fopenmp openmp_demo.cpp -o openmp_demo && ./openmp_demo

在 8 核心 CPU 的计算机上输出：

➜  src_mat git:(master) ./openmp_demo 
Hello, OpenMP!
Hello, OpenMP!
Hello, OpenMP!
Hello, OpenMP!
Hello, OpenMP!
Hello, OpenMP!
Hello, OpenMP!
Hello, OpenMP!
➜  src_mat git:(master)

在开启了 OpenMP 的那一行，一行语句被拆分到 8 个核心分别执行，所以输出了 8 次。

// openmp_demo_for.cpp
#include <cstdio>
#include <omp.h>

int main(int argc, char **argv) {
  #pragma omp parallel
	{
		#pragma omp for
		for (int i = 0; i < 64; i++) printf("i = %d\\n", i);
	}
  return 0;
}

使用 #pragma omp for 对 for 循环开启 OpenMP。 g++ -o openmp_demo_for openmp_demo_for.cpp -fopenmp && ./openmp_demo_for 输出为：

i = 0
i = 1
i = 2
i = 3
i = 4
i = 5
i = 6
i = 7
i = 16
i = 17
i = 18
i = 19
i = 20
i = 21
i = 22
i = 23
i = 8
i = 9
i = 10
i = 11
i = 12
i = 13
i = 14
i = 15
i = 48
......

运行环境为 8 核心 CPU。所以可知，OpenMP 现在按照 8 个线程分组执行。

使用 OpenMP 优化矩阵乘法

代码基本和 lesson04 一样，添加了 OpenMP 操作对比，优化了数据展示方式。

添加的 OpenMP 优化如下：

  1. 普通计算的时候用 OpenMP 优化

    Mat* mat_mul_openmp_native(Mat* a, Mat* b, Mat* c) {
  // 检查是否合法
  if ((!a || !b || !c) || (a->w != b->h)) return NULL;
  // 计时，超时的话就直接返回
  int k = a->w;
#pragma omp parallel
  {
#pragma omp for
    for (int x = 0; x < c->w; x++) {
      for (int y = 0; y < c->h; y++) {
        double sum = 0;
        for (int i = 0; i < k; i++) sum += a->data[x][i] * b->data[i][y];
        c->data[x][y] = sum;
      }
    }
  }
  return c;
}

    如在最内层循环开启 OpenMP，需要使用 #pragma omp parallel for reduction(+:sum) 。

  2. 用 OpenMP 优化 lesson04 中的任务调度程序

    Mat* mat_mul_openmp(Mat* a, Mat* b, Mat* c, int unrolling) {
  // 检查是否合法
  if (a->w != b->h) {
    return NULL;
  }
  // 首先对 b 进行一个置的转
  Mat* t = mat_transpose(b);
  // 初始化任务列表和线程池
  mat_task_tail = a->h * b->w;
  mat_task_list = malloc(sizeof(int*) * mat_task_tail);
  assert(mat_task_list);
  mat_task_data = malloc(sizeof(int) * mat_task_tail * 2);
  assert(mat_task_data);
  // 初始化任务
  for (int x = 0; x < c->h; x++) {
    for (int y = 0; y < c->w; y++) {
      mat_task_list[x * c->w + y] = &mat_task_data[(x * c->w + y) * 2];
      mat_task_list[x * c->w + y][0] = x;
      mat_task_list[x * c->w + y][1] = y;
    }
  }
  // 填满线程参数
  mat_mul_thread_t* thread_data =
      malloc(sizeof(mat_mul_thread_t) * mat_task_tail);
  for (int i = 0; i < mat_task_tail; i++) {
    thread_data[i].a = a;
    thread_data[i].b = t;
    thread_data[i].c = c;
    thread_data[i].id = i;
    thread_data[i].single = 1;
    thread_data[i].single_index = i;
    thread_data[i].unrolling = unrolling;
  }
  #pragma omp parallel for
  for (int t = 0; t < mat_task_tail; t++) {
    mat_mul_cell(thread_data + t);
  }
  return c;
}

    这段代码中首先单个任务的参数生成好，方便调用的时候独立传参。

    lesson04 中是首先填满 CPU 核心大小的线程池，然后每个线程独立领取任务，这里是每个线程只执行一次，不需要领取任务，线程执行之后退出再由 OpenMP 调用执行下一个线程。

实验结果

运行环境标在图片标题上。（CPU 所标频率为默频，运行时频率更高；SIMD 指 AVX256 指令集优化。）

小矩阵情况：$N \in [2^1, \lfloor 2^{7.5} \rfloor]$，重复 50 次取平均值

图1，openmp_update_s1_m7_r50_cubic.png

图1，openmp_update_s1_m7_r50_cubic.png