[Boost-users] OpenMP and Boost - can they be used together?

Subject: [Boost-users] OpenMP and Boost - can they be used together?
From: Max S. Kaznady (max.kaznady_at_[hidden])
Date: 2010-09-15 19:34:57

Hello,

I wrote a simple parallel program which computes the sum of elements
across a matrix using Boost. When I use OpenMP, the program runs give
different sums. When I use just simple arrays like double array[N][N],
instead of boost::numeric::ublas::matrix<double>, the code work just
fine and as expected. Does anyone know what is going on? It looks like
the score function is not being modified properly, because Boost array
returns just fine.

The function compute_score() creates a matrix and computes the sum of
its elements in parallel, copying the result to matrix temp and then
printing temp, to make sure that the results can be copied in parallel
too. The function is run infinitely many times, and if its return
values differ then the program exits.

[code]
#include <iostream>
#include <omp.h>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>

#define ENABLE_OPENMP_SCORE 1
#define N 10
#define NUMTHREADS 2

using std::cout;
using std::endl;
using boost::numeric::ublas::matrix;

double compute_score() {
  matrix<double> data1(N, N);
  for (unsigned i = 0; i < N; ++i)
    for (unsigned j = 0; j < N; ++j)
      data1(i, j) = i * N + j;

  double sc = 0.0;
  unsigned i, j;
  matrix<double> temp(N, N);

#pragma omp parallel \
  shared(sc,temp) \
  private(i,j) \
  if (ENABLE_OPENMP_SCORE)
  {
    omp_set_num_threads(NUMTHREADS);
#pragma omp master
    sc = 0.0;

#pragma omp for reduction(+:sc) schedule(static) collapse(2) nowait
    for (i = 0; i < N; ++i) {
      for (j = 0; j < N; ++j) {
        temp(i, j) = data1(i, j);
        sc += data1(i, j);
      }
    }
  }
  cout << temp << endl;
  return sc;
}

int main(void) {

  double oldval = compute_score(), newval;
  cout.precision(32);
  cout << "--" << endl;
  while (1) {
    newval = compute_score();
    cout << "--" << endl;
    if (newval != oldval) {
      cout << " newval - oldval = diff, " << newval << " - " << oldval << " = "
          << newval - oldval << endl;
      exit(1);
    }
    oldval = newval;
  }

  return 0;
}
[/code]

This problem is driving me insane... the code is really simple, but
for some reason after using Boost it does not function as expected.

Cheers,
Max