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;

  #pragma omp parallel for num_threads(nThreads) reduction(+: s)
  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;

  #pragma omp parallel for num_threads(nThreads) reduction(+: s)
  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()
{
  const short maxThreads = 10;
  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.