#include <stdio.h>
#include <math.h>

#define  n     30
#define  work  10000

void init(double a[n][n]) {
  /* initial values for a: a(i,j) = i+j */

  int  i,j;
    
  for (i=0; i<n; i++) {
    for (j=0; j<n; j++) {
      a[i][j] = i + j + 2;
    }
  }
}

void check(double v[n]) {
  /* check results */
  /* here: print total sum */
    
  double  sum1, sum2;
  int     i;
    
  sum1 = 0;
  for (i=0; i<n; i++) {
    sum1 += v[i];
  }

  sum2 = n*n*n + n*n;
    
  printf("sum is       : %lf\n", sum1);
  printf("and should be: %lf\n", sum2);
}


double f(double x, int w) {
  /* time consuming way of setting f = x + eps */

  int    i;
  double t1, t2;
    
  t1 = x;
  for (i=0; i<w; i++) {
    t1 = cos(t1);
  }
    
  t2 = x;
  for (i=0; i<w; i++) {
    t2 = sin(t2);
  }
  
  return x + t1 + t2 - 0.7390851332151607;
}


int main(int argc, char *argv[]) {

  double a[n][n];
  double v[n];
  int    i, j;

  for (i=0; i<n; i++) {
    v[i] = 0.0;
  }
  init(a);   /* initial values for a */

#pragma omp parallel
{
#pragma omp for
  for (i=0; i<n; i++) {
    for (j=0; j<n; j++) {
      v[i] = v[i] + f(a[i][j], work);
    }
  }
}

  check(v);  /* check result */
}