使用openmp进行共享内存编程

预处理指令pragma:在系统中加入预处理器指令一般是用来允许不是基本c语言规范部分的行为。不支持pragma的编译器会忽略pragma指令提示的那些语句，这样就允许使用pragma的程序在不支持它们的平台上运行。

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

void Hello(void); // Thread function

int main(int argc, char* argv[]) {
  // Get number of threads from command line
  int thread_count = strtol(argv[1], NULL, 10);
#pragma omp parallel num_threads(thread_count)
  Hello();

  return 0;
}

void Hello(void) {
int my_rank = omp_get_thread_num();
int thread_count = omp_get_num_threads();

printf("Hello from thread %d of %dnn
", my_rank, thread_count);

} //Hello

例子的分析：

# pragma omp parallel  :

最基本的并行原语，用于运行代码块的线程数可以动态生成。

当程序到达parallel指令时，原来的线程继续执行，另外的线程被启动。在openmp语法中，执行并行块的线程集合（原始线程和新的线程被称为线程组，原始的线程被称为主线程，额外的线程称为从线程。每个线程组成员都调用指令后的代码块。

# pragma omp parallel num_threads ( thread_count )

num_threads（）：一个从句例子（用于修饰原语），可用于指定线程数量

#include <omp.h>

使用openmp必须含omp.h头文件

long strtol(const char* number p,char** end p,int base)

使用stdlib.h中的strtol来获得线程数

ps:一些系统因素可能会限制可以启动的线程数量；OpenMP 并不保证能够启动指定个的线程；

多数系统能够启动上百甚至上千的线程；除非启动的太多，一般都能满足要求。

eg:梯形积分法

如果每个子区间有相同的宽度，并且定义h=(b-a)/n，xi=a+ih，i=0, 1, ...， n,那么近似值将是：

 //串行算法实现
// Input: a, b, n 
h = (b*a)/n;
approx = (f(a) + f(b))/2.0;
for (i = 1; i <= n-1; i++) {
x_i = a + ih;
approx += f(x_i);
}
approx = h*approx;

第一种尝试：1) 定义两种类型的任务:

　　　　a) 计算单个梯形的面积；

　　　　b) 将面积加起来。

2) 在第一阶段，没有通信开销；但第二阶段每个任务需要通信。考虑一个问题：结果不可预估——引入互斥量

# pragma omp critical
   global_result += my_result ;

第一个版本：

 #include <stdio.h>
 #include <stdlib.h>
 #include <omp.h>

  void Trap(double a, double b, int n, double  global_result p);

  int main(int argc, char  argv[]){
    double  global_result = 0.0;
    double  a, b;
    int     n;
    int     thread_count;

    thread_count = strtol(argv[1], NULL, 10);
    printf("Enter a, b, and n n");
    scanf("%lf %lf %d", &a, &b, &n);
#   pragma omp parallel num_threads(thread_count)
    Trap(a, b, n, &global_result);

    printf("With n = %d trapezoids, our estimate n", n);
    printf("of the integral from %f to %f = %.14e n",
    a, b, global_result);
    return 0;
}    /∗  main ∗/
    
    void Trap(double a, double b, int n, double* global_result_p)
      double  h, x, my_result;
      double  local_a, local_b;
      int  i, local n;
      int my_rank = omp_get_thread_num();
      int thread_count = omp_get_num_threads();
        
      h = (b−a)/n;
      local_n = n/thread_count;
      local_a = a + my_rank*local_n*h;
      local_b = local_a + local_n*h;
      my_result = (f(local_a) + f(local_b))/2.0;
      for (i = 1; i <= local_n−1; i++){
        x = local_a + i*h;
        my_result += f(x);
       }
  `   ` my_result = my_result*h;
   #  pragma omp critical
      ∗global_result_p += my_result;
  }    /∗  Trap ∗/

规约从句：

替代（在parallel块中声明一个私有变量和将临界区移到函数调用之）

归约：将相同的归约操作符重复的应用到操作数序列来得到一个结果的计算。

所有操作的中间结果存储在一个变量中：归约变量

reduction(<operator>:<variable list>)

新的代码：

global_result = 0.0;
#  pragma omp parallel num threads(thread count)
 reduction(+: global_result)
global_result += Local_trap(double a, double b, int n);

parallel for：

能够生成一队线程来执行接下来的语句块；

语句块必须是一个for循环；

通过将循环切分给不同的线程来实现并行。

只有迭代次数确定的循环才可以被并行化。

h = (b−a)/n;
approx = (f(a) + f(b))/2.0;
#  pragma omp parallel for num threads(thread_count) reduction(+: approx)
for (i = 1; i <= n−1; i++)
approx += f(a + i∗h); approx = h∗approx;

可被并行化的for循环形式：

ps: index 必须是整数或者指针 (e.g., 不能是浮点数)；

start, end, 和 incr 必须具有相应的类型。例如, 如果index 是一个指针, 那么 incr 必须是一个整型；

start, end, 和 incr 在循环执行过程中不能被修改；

在循环执行过程中, 变量 index 只能被for语句修改。

数据依赖：1.OpenMP 编译器并不检查循环迭代中的数据依赖问题；

2.一般来说，OpenMP无法处理带有数据依赖的循环。

解决思路：设计私有变量并且保证其私有作用域（private子句）

default子句：编译器强制要求程序员指定在块中使用的外部变量的作用范围。

double sum = 0.0;
# pragma omp parallel for num threads(thread count)
default(none) reduction(+:sum) private(k, factor)
 shared(n)
for (k = 0; k < n; k++){
  if (k % 2 == 0)
    factor = 1.0;
  else
    factor = −1.0;
  sum += factor/(2∗k+1);
}

for指令：并不创建线程，使用已经在parallel块中创建的线程。

#  pragma omp for

解决循环调用问题：

schedule ( type , chunksize )
//chunksize 是一个正整数

type 可以是:

static: 提前把任务分配好；

dynamic or guided: 在运行时动态分配；

dynamic:任务被分成 chunksize 大小的连续段；每个线程执行一小块, 当有一个线程执行完时, 它会请求获得1个新的；重复上述过程，直到完成计算；

chunksize 可以被去掉;当去掉时, chunksize 默认为1.

guided:每个线程执行一小块, 当有一个线程执行完时, 它会请求获得1个新的；但是，新的任务块是不断变小的；如果不指定chunksize，那么默认会降到1.如果指定了chunksize, 则会降到指定的chunksize, 除了最后一块可能小于chunksize.

auto: 编译器或者运行时系统决定调度策略；

runtime: 运行时决定。