Хочу распараллелить алгоритм нахождения интеграла с помощью tbb. Но похоже я что-то неправильно делаю. Функция IntegralTBB возвращает каждый раз новое значение, похоже гонка данных. Как исправить?
#include <iostream>
#include "tbb/tbb.h"
double IntegralSeq(const double a1, const double b1,
const double a2, const double b2,
const double h)
{
int n1, n2;
// локальная переменная
// для подсчета интеграла
double sum = 0.0f;
// координата точки сетки по оси x
// координата точки сетки по оси y
double x, y;
// количество точек сетки интегрирования
// n1 - по координате x
// n2 - по координате y
n1 = (int)(b1 - a1) / h;
n2 = (int)(b2 - a2) / h;
for(int i = 0; i < n1; ++i)
{
for(int j = 0; j < n2; ++j)
{
//вычисление координат точки
x = a1 + i * h + h / 2;
y = a2 + j * h + h / 2;
//вычисление интеграла
sum += ((exp(sin(x * M_PI) * cos(y * M_PI)) + 1) /
((b1 - a1) * (b2 - a2))) * h * h;
}
}
return sum;
}
double IntegralTBB(const double a1, const double b1,
const double a2, const double b2,
const double h)
{
int n1, n2;
n1 = (int)(b1 - a1) / h;
n2 = (int)(b2 - a2) / h;
return tbb::parallel_reduce(tbb::blocked_range2d<int>(0,n1,0,n2),0.f,[&](const tbb::blocked_range2d<int> & r, float sum)
{
for(int i = r.rows().begin(); i != r.rows().end(); ++i) {
for(int j = r.cols().begin(); j != r.cols().end() ; ++j)
{
double x = a1 + i * h + h / 2;
double y = a2 + j * h + h / 2;
sum += (exp(sin(x * M_PI) * cos(y * M_PI)) + 1) /
((b1 - a1)*(b2 - a2) * h * h);
}
}
return sum;
},std::plus<int>());
}
int main()
{
tbb::tick_count t0 = tbb::tick_count::now();
double res = IntegralSeq(0.0f, 6.0f, 0.0, 6.0, 0.001);
tbb::tick_count t1 = tbb::tick_count::now();
std:: cout << " result seq:" << res << " time:" << (t1-t0).seconds();
t0 = tbb::tick_count::now();
res = IntegralTBB(0.0f, 6.0f, 0.0, 6.0, 0.001);
t1 = tbb::tick_count::now();
std:: cout << " result tbb:" << res << " time:" << (t1-t0).seconds();
return 0;
}
std::plus<int>()
, а должно бытьfloat
... в остальном вроде всё в порядке.... и раз уж все вычисления проходят в double, то сумму я тоже сделал бы в double... – Fat-Zer 24 фев '18 в 23:47n1 = (int) ((b1 - a1) / h);
? – AnT 25 фев '18 в 0:41h * h
в числителе, а во второй - в знаменателе. – AnT 25 фев '18 в 0:54