# Parallel numerical integration with OpenMP

Numerical integration is the approximate computation of an integral using numerical techniques. Today, I am going to implement two simple algorithms, which uses Rectangle and Trapezoidal Method.

### What is OpenMP?

OpenMP is an API, which supports C, C++ and Fortran. It provides portable and scalable model for developers of shared memory parallel applications. It helps us to focus only on algorithm implementation (in most cases).

### Let’s start coding

So, firstly – we have to create a function, which will return output (`y`) of our input (`x`) for specific formula. I choosed Sine.

``````
double f(const double x)
{
return sin(x);
}
``````

Since we had a function, we are able to prepare Rectangle Method.

``````
const Result rectangleMethod(const double x1, const double x2, const double dx, const int nThreads)
{
const int N = static_cast<int>((x2 - x1) / dx);
double now = omp_get_wtime();
double s = 0;

for (int i = 1; i <= N; i++) s += f(x1 + i * dx);

s *= dx;

return { omp_get_wtime() - now, s };
}
``````

We can see that our code looks like just sequential implementation of algorithm. A lot of stuff is done by OpenMP and we don’t have to worry about it. But still, some things may be a little bit unclear here. Firstly, how `reduction` clause works?

Reduction clause. A private copy for each list variable is created for each thread. At the end of the reduction, the reduction variable is applied to all private copies of the shared variable, and the final result is written to the global shared variable. – documentation

Second thing – how `Result` struct looks?

``````
struct Result
{
double timestamp, area;
};
``````

And that’s it. Now, it is time for second method, which uses Trapezoidal rule.

``````
const Result trapezoidalMethod(const double x1, const double x2, const double dx, const int nThreads)
{
const int N = static_cast<int>((x2 - x1) / dx);
double now = omp_get_wtime();
double s = 0;

for (int i = 1; i < N; i++) s += f(x1 + i * dx);

s = (s + (f(x1) + f(x2)) / 2) * dx;

return { omp_get_wtime() - now, s };
}
``````

Looks very similar, almost same, right? It is incredible, how it is easy to work with OpenMP. The last thing – `main()` function:

``````
int main()
{
short method;
double x1, x2, dx;

cout << fixed << setprecision(8);
try
{
while (true)
{
cout << "X1: "; cin >> x1;
cout << "X2: "; cin >> x2;
cout << "dx: "; cin >> dx;
cout << "Method (1 - rectangle, 2 - trapezoidal): "; cin >> method;

list<pair<short, Result>> results;
for (int i = 0; i < maxThreads; i++)
{
Result result = (method == 1) ?
rectangleMethod(x1, x2, dx, i + 1) :
trapezoidalMethod(x1, x2, dx, i + 1);

pair<short, Result> s_result(i + 1, result);
results.push_back(s_result);
}

cout << endl << "Results:" << endl;
for (auto & result : results)
{
cout << "Threads: " << result.first;
cout << ", timestamp: " << result.second.timestamp;
cout << ", area: " << result.second.area << endl;
}
cout << endl;
}
}
catch (exception & e)
{
cout << e.what() << endl;
}
cin.get();
return 0;
}
``````

Everything is done. Finally, take a look how our solution works.

Full source code available at GitHub.