OpenMP学习 第七章 OpenMP任务
第七章 OpenMP任务
不规则问题
总结前面的内容,我们所学内容大都符合下面的逻辑:在串行程序中找到计算密集型循环,添加并行共享工作循环构造将其转化为并行应用程序.
但是,实际上,有一类重要的问题并不规则.它并不直接映射到嵌套循环索引上;亦或者即使能够映射,每个循环迭代的工作也是千变万化的.这类问题统称为 不规则问题.
不规则问题 很多,它们包括基于稀疏数据结构的问题,其工作是高度可变的,或者是控制流和依赖性不可预测的问题,这些问题非常复杂,无法用规则迭代空间来表示.
为了解决不规则问题的并行化,OpenMP通过 任务(tasks) 来实现目标.
在任务之前,我们组织语言的中心实体是线程.有了任务,我们现在就有了独立于线程的工作单元.
任务
single构造 表示一个只由一个线程执行的代码块,该线程为线程组中的任何线程.
- 创建single共享工作构造的语句:
#pragma omp single[nowait]
{
//the body of single
}
下面给出一个使用single构造的例子,这是在处理OpenMP/MPI混合程序时的常见模式:由于MPI的某些实现不太好与多个线程合作,因而可以指定一个线程与MPI合作,以在不同迭代MPI进程间交换边界.
#pragma omp parallel
{
do_many_things();
#pragma omp single
{
exchange_boundaries();
}
do_many_other_things();
}
- 创造一个显式任务的构造:
#pragma omp task[clause[,[clause]]...]
{
//the body of task
}
- 强制等待前文显式任务:
#pragma omp taskwait
当一个线程遇到一个任务构造时,它可能会立即执行任务,也可能会将其推迟到以后执行.正是这种延迟执行使得任务变得有趣.
类似于barrier,我们也可以通过taskwait显式声明等待前文所提及任务的执行.
为了帮助理解,我们通过下面这个名为: 薛定谔程序 的程序来加深理解.
#include <iostream>
#include <vector>
#include <cstdlib>
#include <ctime>
#include <windows.h>
#include <omp.h>
import <format>;
#define TURNS 10
#define MAX_DELAY_TIMES 1000
int main()
{
std::vector<double>times;
double temp_time, run_time;
int count = 0;
srand(time(0));
omp_set_num_threads(2);
count = 0;
for (int iter = 0; iter < TURNS; iter++) {
bool is_alive;
time_t rand_1{ rand() % MAX_DELAY_TIMES };
time_t rand_2{ rand() % MAX_DELAY_TIMES };
#pragma omp parallel barrier
temp_time = omp_get_wtime();
#pragma omp parallel
{
#pragma omp task
{
Sleep(rand_1);
is_alive = false;
}
#pragma omp task
{
Sleep(rand_2);
is_alive = true;
}
}
if (is_alive)
count++;
run_time = omp_get_wtime() - temp_time;
times.push_back(run_time);
}
double aver{ 0.0 };
for (auto iter : times)
aver += iter;
aver /= TURNS;
std::cout << std::format("the average time is: {:.15f},the count is: {}", aver, count);
return 0;
}
理论上来说,该程序应当是有随机结果的(如果windows.h中的Sleep能够使得OpenMP中单个thread休眠).
在处理is_alive处,这里呈现了所谓的"良性数据竞争".
显式任务构造非常灵活,可以以多种方式使用.通常的一种使用模式为: 选择一个线程通过single构造来创建任务,而组中其他线程在栅栏处等待.
#pragma omp parallel
{
#pragma omp single
{
#pragma omp task
task_1();
#pragma omp task
task_2();
#pragma omp task
task_3();
}//end of single region
}//end of parallel region
任务是不规则问题的理想选择.一个面向任务的系统被称为 自动平衡负载.
任务的数据环境
任务的数据环境规则与其他OpenMP构造的规则类似,有两个区别:
- 数据环境绑定的是任务,而不是遇到任务的线程.
- 如果一个变量在遇到任务时是私有的,那么它将被默认为firstprivate.
即可以进一步解释为下面四条:
- 数据环境绑定的是任务,而不是遇到任务的线程.
- 如果一个变量在任务构造上是shared,那么在构造内部其为shared
- 如果一个变量在任务构造上是private,那么在构造内部其为firstprivate
- 如果一个变量在任务构造上是firstprivate,那么在构造内部其为firstprivate
利用任务的基础设计模式
- 链表遍历
为了探究使用任务来解决并行问题的设计模式,我们从一个简单的链表程序开始:
问题描述: 下面是一段代码,processWork函数代表一个工作块,代码的功能为遍历一个链表.
p=head;
while(p!=NULL){
processWork(p);
p=p->next;
}
为了解决该问题,首先思考该问题的串行解法:
#include <iostream>
#include <cstdlib>
#include <omp.h>
#define NODE_NUM 10
typedef struct node {
int data;
int procResult;
struct node* next;
node():data(0), procResult(0), next(nullptr){}
}Node,* List;
void initList(List p)
{
Node* root{ p };
Node* temp_node;
p->data = 0;
for (int i = 1; i < NODE_NUM; i++) {
temp_node = new Node;
temp_node->data = i;
root->next = temp_node;
root = temp_node;
}
return;
}
void processWork(Node* n)
{
n->procResult = (n->data * n->data);
return;
}
void deleteList(List p)
{
Node* temp_node = p->next;
for (; p != temp_node;) {
temp_node = p;
while (temp_node->next != nullptr && temp_node->next->next != nullptr)
temp_node = temp_node->next;
delete temp_node->next;
temp_node->next = nullptr;
}
delete p;
return;
}
int main()
{
List list = new Node;
initList(list);
int count = 0;
Node* head = list;
while (head != nullptr) {
processWork(head);
count++;
head = head->next;
}
deleteList(list);
return 0;
}
为了让其并行化,我们首先尝试不使用任务,而是用共享任务循环构造解决该问题.这需要借助一个数组.
#include <iostream>
#include <cstdlib>
#include <omp.h>
#define NODE_NUM 10
#define CHUNK 2
typedef struct node {
int data;
int procResult;
struct node* next;
node():data(0), procResult(0), next(nullptr){}
}Node,* List;
void initList(List p)
{
Node* root{ p };
Node* temp_node;
p->data = 0;
for (int i = 1; i < NODE_NUM; i++) {
temp_node = new Node;
temp_node->data = i;
root->next = temp_node;
root = temp_node;
}
return;
}
void processWork(Node* n)
{
n->procResult = (n->data * n->data);
return;
}
void deleteList(List p)
{
Node* temp_node = p->next;
for (; p != temp_node;) {
temp_node = p;
while (temp_node->next != nullptr && temp_node->next->next != nullptr)
temp_node = temp_node->next;
delete temp_node->next;
temp_node->next = nullptr;
}
delete p;
return;
}
int main()
{
List list = new Node;
Node** parr = new Node*[NODE_NUM];
initList(list);
int count = 0;
Node* head = list;
while (head != nullptr) {
parr[count] = head;
head = head->next;
count++;
}
#pragma omp parallel for schedule(static, CHUNK)
for (int i = 0; i < NODE_NUM; i++)
processWork(parr[i]);
deleteList(list);
return 0;
}
而后,为了方便比较,我们再使用single任务构造对其进行并行化.
#include <iostream>
#include <cstdlib>
#include <omp.h>
#define NODE_NUM 10
#define CHUNK 2
typedef struct node {
int data;
int procResult;
struct node* next;
node():data(0), procResult(0), next(nullptr){}
}Node,* List;
void initList(List p)
{
Node* root{ p };
Node* temp_node;
p->data = 0;
for (int i = 1; i < NODE_NUM; i++) {
temp_node = new Node;
temp_node->data = i;
root->next = temp_node;
root = temp_node;
}
return;
}
void processWork(Node* n)
{
n->procResult = (n->data * n->data);
return;
}
void deleteList(List p)
{
Node* temp_node = p->next;
for (; p != temp_node;) {
temp_node = p;
while (temp_node->next != nullptr && temp_node->next->next != nullptr)
temp_node = temp_node->next;
delete temp_node->next;
temp_node->next = nullptr;
}
delete p;
return;
}
int main()
{
List list = new Node;
Node** parr = new Node*[NODE_NUM];
initList(list);
Node* p;
#pragma omp parallel
{
#pragma omp single
{
p = list;
while (p != nullptr)
{
#pragma omp task firstprivate(p)
{
processWork(p);
}//end of task creation
p = p->next;
}
}//end of single region
}//end of parallel region
deleteList(list);
return 0;
}
可以看到,通过显式任务single构造,我们十分优雅地解决了链表遍历的并行化.
在这个过程中,单个线程首先抓取列表的头部,然后用while循环遍历链表,并将每个节点的进程打包成一个任务,然后继续下一个结点,直至到达列表尾部.
- Fibonacci问题
斐波那契问题是算法问题中的一个经典问题.在并行程序设计领域中更是"Hello World"级别程序的存在
问题描述: 斐波那契数列被定义为:
$$F(n)=F(n-1)+F(n-2)$$
同时有:$F(0)=0,F(1)=1$.请对该问题进行并行化.
为了解决该问题,首先做出串行算法的解答:
#define N 100
int fib(int n)
{
if(n<2)
return n;
return fib(n-1)+fib(n-2);
}
int main()
{
fib(N);
return 0;
}
这是一个$O(n^2)$的递归实现,其效率很低,但是却符合并行的思路.
#define N 100
int fib(int n)
{
int x,y;
if(n<2)
return n;
#pragma omp task shared(x)
x=fib(n-1);
#pragma omp task shared(y)
y=fib(n-2);
#pragma omp taskwait
return x+y;
}
int main()
{
#pragma omp parallel
{
#pragma omp single
fib(N);
}//end of parallel region
return 0;
}
这个问题的解法非常经典,因为它很好地体现了并行算法分而治之的思想.
分而治之模式 是递归程序的基本模式,就像循环级并行模式是基于网格规则网络代码的基本模式一样.
在分而治之模式向,我们递归地将问题分割为越来越小的子问题,直到子问题变得足够小,以至于直接解决他们就可以.
因而,在实现分而治之算法时,需要做的一个重要的决定就是何时进行直接求解.
- Pi循环问题:
问题描述: 数学上,我们已经知道了:
$$\int\limits1_0\dfrac{4.0}{(1+x2)}dx=\pi$$
为了对其进行计算,将积分近似为多个矩形面积的和:
$$\sum\limits^N_{i=0}F(x_i)\Delta x\approx\pi$$
于是,现在需要设计一个并行算法对其进行解决.
在这里,我们回顾我们先前在第四章所学的利用中点规则进行数值积分的解法https://www.cnblogs.com/mesonoxian/p/17968788
#include <iostream>
#include <omp.h>
import <format>;
static long num_steps = 100000000;
double step;
int main()
{
double x, pi, sum = 0.0;
double start_time, run_time;
step = 1.0 / (double)num_steps;
start_time = omp_get_wtime();
for (int i = 0; i < num_steps; i++) {
x = (i + 0.5) * step;
sum += 4.0 / (1.0 + x * x);
}
pi = step * sum;
run_time = omp_get_wtime() - start_time;
std::cout << std::format("pi={},{} steps,{} secs\n",
pi,
num_steps,
run_time
);
}
现在我们希望将其能够通过 分而治之模式 求解.因而对我们来说,重要的部分还有两步:
- 将问题分解为子问题
- 决定何时进行直接求解
为了进行递归求解,首先将问题求解的主要部分代码单独取出为一个函数:
double pi_cal(int Nstart, int Nfinish, double step)
{
double x, sum = 0.0;
//#pragma omp for reduction(+:sum)
for (int i = Nstart; i < Nfinish; i++) {
x = (i + 0.5) * step;
sum += 4.0 / (1.0 + x * x);
}
return sum;
}
然后将问题分解为子问题,思考如何将其递归化:
double pi_recur(int Nstart, int Nfinish, double step)
{
long blength = Nfinish - Nstart;
double sum_1, sum_2;
#pragma omp task shared(sum_1)
sum_1 = pi_comp(Nstart, Nfinish - blength / 2, step);
#pragma omp task shared(sum_2)
sum_2 = pi_comp(Nfinish - blength / 2, Nfinish, step);
#pragma omp taskwait
return sum_1 + sum_2;
}
将其组合起来,设置递归结束条件(即何时进行直接运算):
double pi_comp(int Nstart,int Nfinish,double step)
{
double sum;
if(Nfinish-Nstart<MIN_BLK)
sum = pi_cal(step);
else
sum = pi_recur(Nstart, Nfinish, step);
return sum;
}
将其归于原位,于是便顺利将其分解为递归形式的子问题:
再通过观察易发现,在pi_cal处可以进行并行循环共享工作构造的归约,在pi_recur处可以进行分而治之模式的并行.
将两处并行,再对程序框架进行微调,得到了最终并行化的结果:
#include <iostream>
#include <fstream>
#include <omp.h>
import <format>;
#define TURNS 100
#define MIN_BLK 1024 * 256
#define PI 3.141592653589793
long num_steps = 1024 * 1024 * 1024;
double step;
double pi_cal(int Nstart, int Nfinish, double step);
double pi_recur(int Nstart, int Nfinish, double step);
double pi_comp(int Nstart, int Nfinish, double step);
int main()
{
std::ofstream out;
out.open("example.csv", std::ios::ate);
out << "NTHREADS,pi,err,run_time,num_steps" << std::endl;
double sum = 0.0;
for (int NTHREADS = 1; NTHREADS < TURNS; NTHREADS++) {
double start_time, run_time;
double pi, err;
pi = sum = 0.0;
int actual_nthreads;
step = 1.0 / (double)num_steps;
omp_set_num_threads(NTHREADS);
start_time = omp_get_wtime();
actual_nthreads = omp_get_num_threads();
#pragma omp parallel
{
#pragma omp parallel single
sum = pi_comp(0, num_steps, step);
}//end of parallel
pi = step * sum;
err = pi - PI;
run_time = omp_get_wtime() - start_time;
std::cout << std::format("pi is {} in {} seconds {} thrds.step is {},err is {}",
pi,
run_time,
actual_nthreads,
step,
err
) << std::endl;
out << std::format("{},{:.15f},{:.15f},{:.15f},{}",
NTHREADS,
pi,
err,
run_time,
num_steps
) << std::endl;
}
out.close();
return 0;
}
double pi_cal(int Nstart, int Nfinish, double step)
{
double x, sum = 0.0;
#pragma omp for reduction(+:sum)
for (int i = Nstart; i < Nfinish; i++) {
x = (i + 0.5) * step;
sum += 4.0 / (1.0 + x * x);
}
return sum;
}
double pi_recur(int Nstart, int Nfinish, double step)
{
long blength = Nfinish - Nstart;
double sum_1, sum_2;
#pragma omp task shared(sum_1)
sum_1 = pi_comp(Nstart, Nfinish - blength / 2, step);
#pragma omp task shared(sum_2)
sum_2 = pi_comp(Nfinish - blength / 2, Nfinish, step);
#pragma omp taskwait
return sum_1 + sum_2;
}
double pi_comp(int Nstart, int Nfinish, double step)
{
double sum;
if (Nfinish - Nstart < MIN_BLK)
sum = pi_cal(Nstart, Nfinish, step);
else
sum = pi_recur(Nstart, Nfinish, step);
return sum;
}
经过实验,我们发现如此得到的结果精确度已经达到了非常高的水平,在double数据类型表示范围内误差为0,且运行的时间也在我们接受范围内.
