簡(jiǎn)述

會(huì)用多個(gè)版本來寫。操作是用命令行來操作。
除了第二版本是在內(nèi)容上基于第一個(gè)版本的完善，其他都是降低算法的復(fù)雜程度的，可以放心閱讀

文章目錄

• • 簡(jiǎn)述
  • 版本1
  • 版本2
  • 版本3
  • 版本4

版本1

• 要求n被thread_count整除 不然可以會(huì)有計(jì)算錯(cuò)誤。

#include <iostream> #include <omp.h>#define FUN(x) (x * x)using namespace std; #pragma warning(disable : 4996) void Trap(double a, double b, int n, double *global_result_p); int main(int argc, char **argv) {if (argc < 5) return 0;double global_result = 0.0;double a, b;int n;int thread_count = strtol(argv[1], NULL, 10);a = strtod(argv[2], NULL);b = strtod(argv[3], NULL);n = strtol(argv[4], NULL, 10);// n 必須被 thread_count 整除 #pragma omp parallel num_threads (thread_count)Trap(a, b, n, &global_result);cout << "With n = " << n << " trapzoids, our estimate\n";cout << "of the integral from " << a << " to " << b <<" = "<< global_result<< endl; }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; // same in all threadslocal_a = a + my_rank * local_n * h;local_b = local_a + local_n * h;x = local_a; // initial xmy_result = (FUN(local_a) + FUN(local_b)) / 2;for (i = 1; i < local_n; ++i) {x += h;my_result += FUN(x);}my_result *= h; # pragma omp critical *global_result_p += my_result; }

編譯好源文件之后，就用下面的命令執(zhí)行

• 這里是用4個(gè)線程，來計(jì)算在0，1區(qū)間上的$x^2$的用n=64劃分的梯形積分公式。結(jié)果接近1/3 可以通過提高n來提高精度。但這個(gè)版本只能計(jì)算n為thread_count的整數(shù)倍的時(shí)候

PS D:\C++\VS\repo\OpenMP-TEST\Debug> ./OpenMP-TEST 4 0 1 64 With n = 64 trapzoids, our estimate of the integral from 0 to 1 = 0.333374

版本2

• 可以使用不被thread_count整除的n
• 改進(jìn)代碼： 添加了一個(gè)變量left表示剩余。剩余的讓前l(fā)eft個(gè)線程每個(gè)再多算點(diǎn)。（一般來說線程數(shù)目都會(huì)比n要小很多。（不然在算法上就會(huì)不精確（n不夠大的話）。））

left = n % thread_count; if (my_rank < left) local_n += 1; if (my_rank < left){ local_a = a + my_rank * local_n * h; } else { local_a = a + left * (local_n + 1) * h + (my_rank - left) * local_n * h; }

實(shí)際的代碼:

#include <iostream> #include <omp.h>#define FUN(x) (x * x)using namespace std; #pragma warning(disable : 4996) void Trap(double a, double b, int n, double *global_result_p); int main(int argc, char **argv) {if (argc < 5) return 0;double global_result = 0.0;double a, b;int n;int thread_count = strtol(argv[1], NULL, 10);a = strtod(argv[2], NULL);b = strtod(argv[3], NULL);n = strtol(argv[4], NULL, 10);#pragma omp parallel num_threads (thread_count)Trap(a, b, n, &global_result);cout << "With n = " << n << " trapzoids, our estimate\n";cout << "of the integral from " << a << " to " << b <<" = "<< global_result<< endl; }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, left;int my_rank = omp_get_thread_num();int thread_count = omp_get_num_threads();h = (b - a) / n;local_n = n / thread_count; // same in all threadsleft = n % thread_count;if (my_rank < left) local_n += 1;if (my_rank < left){ local_a = a + my_rank * local_n * h; }else { local_a = a + left * (local_n + 1) * h + (my_rank - left) * local_n * h; }local_b = local_a + local_n * h;x = local_a; // initial xmy_result = (FUN(local_a) + FUN(local_b)) / 2;for (i = 1; i < local_n; ++i) {x += h;my_result += FUN(x);}my_result *= h; # pragma omp critical *global_result_p += my_result; }

發(fā)現(xiàn)精度不斷的提升，在只改變n的數(shù)目的情況下。

PS D:\C++\VS\repo\OpenMP-TEST\Debug> ./OpenMP-TEST 4 0 1 64 With n = 64 trapzoids, our estimate of the integral from 0 to 1 = 0.333374 PS D:\C++\VS\repo\OpenMP-TEST\Debug> ./OpenMP-TEST 4 0 1 65 With n = 65 trapzoids, our estimate of the integral from 0 to 1 = 0.333373 PS D:\C++\VS\repo\OpenMP-TEST\Debug> ./OpenMP-TEST 4 0 1 66 With n = 66 trapzoids, our estimate of the integral from 0 to 1 = 0.333372 PS D:\C++\VS\repo\OpenMP-TEST\Debug> ./OpenMP-TEST 4 0 1 67 With n = 67 trapzoids, our estimate of the integral from 0 to 1 = 0.33337

版本3

• 對(duì)于部分對(duì)指針不太熟的人來說，我們可以做出下面的改進(jìn)
• 修改部分：

#pragma omp parallel num_threads (thread_count){double my_result = Trap(a, b, n); #pragma omp critical global_result += my_result;} #include <iostream> #include <omp.h>#define FUN(x) (x * x)using namespace std; #pragma warning(disable : 4996) double Trap(double a, double b, int n); int main(int argc, char **argv) {if (argc < 5) return 0;double global_result = 0.0;double a, b;int n;int thread_count = strtol(argv[1], NULL, 10);a = strtod(argv[2], NULL);b = strtod(argv[3], NULL);n = strtol(argv[4], NULL, 10);#pragma omp parallel num_threads (thread_count){double my_result = Trap(a, b, n); #pragma omp critical global_result += my_result;}cout << "With n = " << n << " trapzoids, our estimate\n";cout << "of the integral from " << a << " to " << b <<" = "<< global_result<< endl; }double Trap(double a, double b, int n) {double h, x, my_result;double local_a, local_b;int i, local_n, left;int my_rank = omp_get_thread_num();int thread_count = omp_get_num_threads();h = (b - a) / n;local_n = n / thread_count; // same in all threadsleft = n % thread_count;if (my_rank < left) local_n += 1;if (my_rank < left){ local_a = a + my_rank * local_n * h; }else { local_a = a + left * (local_n + 1) * h + (my_rank - left) * local_n * h; }local_b = local_a + local_n * h;x = local_a; // initial xmy_result = (FUN(local_a) + FUN(local_b)) / 2;for (i = 1; i < local_n; ++i) {x += h;my_result += FUN(x);}my_result *= h;return my_result; }

結(jié)果

PS D:\C++\VS\repo\OpenMP-TEST\Debug> ./OpenMP-TEST 4 0 1 67 With n = 67 trapzoids, our estimate of the integral from 0 to 1 = 0.33337 PS D:\C++\VS\repo\OpenMP-TEST\Debug> ./OpenMP-TEST 4 0 1 68 With n = 68 trapzoids, our estimate of the integral from 0 to 1 = 0.333369 PS D:\C++\VS\repo\OpenMP-TEST\Debug>

版本4

• 使用規(guī)約變量和規(guī)約操作
• 會(huì)讓代碼更加簡(jiǎn)單

#pragma omp parallel num_threads (thread_count) reduction(+:global_result){global_result += Trap(a, b, n);}

• 跟版本3是一樣的。就是讓代碼更加簡(jiǎn)單了而已，用了map-reduce的思路
• reduction(<operator>: <variable list>) 前面的放操作符，后面的設(shè)置變量。
• 但是有需要注意的，如果是減法就需要小心了。在每個(gè)線程上確實(shí)是使用了減法，但是全局上也是可能做了減法，(例如 -1 - (-5)= 4）這樣的事情是有可能發(fā)生的。一般來說，滿足交換律的都不用擔(dān)心這點(diǎn)hh

《新程序員》：云原生和全面數(shù)字化實(shí)踐50位技術(shù)專家共同創(chuàng)作，文字、視頻、音頻交互閱讀

總結(jié)

以上是生活随笔為你收集整理的梯形积分法【OpenMP实现】多个版本的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。