title: Parallel computing for prefix sum using openMP (c code)
body: <p>Have some problems with assigning parallel algorithm to prefix sum issue. I am using openMP for parallel implementation. I have the code in c as below. </p>
  
  <p>Result showing:</p>
  
  <pre><code>seqsum[6] = 28 != parallelsum[6] = 34
  </code></pre>
  
  <p>Please advise. Thanks.</p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;time.h&gt;
  #include "omp.h"
  #include &lt;string.h&gt;
  #define N 10 //33554432 // 2 ^ 25
  #define NUM_THREADS 4 
  
  void computeparallelprefix(int *iplist, int *_pprefixsum, unsigned long size)
  {
    int nthr, *z, *x = _pprefixsum;
    int i, j, tid, work, lo, hi;
  #pragma omp parallel shared(nthr,x,z) private(i,j,tid,work,lo,hi)
    {
      int prev_sum;
      memcpy((void *)x, (void *)iplist, sizeof(int)*size);
  
      // Assume nthr = 2^k
  #pragma omp single
      {
        nthr = omp_get_num_threads();
        z = malloc(sizeof(int)*nthr);
      }   
      tid = omp_get_thread_num();
      work = size /nthr + (i = tid &lt; size%nthr ? 1 : 0);
      lo = (size/nthr)*tid + (i==1 ? tid : size%nthr);
      hi = lo + work;
      if (hi &gt; size)
        hi = size;
  
      // local prefix sum over x
      for(i=lo+1; i&lt;hi; i++)
        x[i] += x[i-1];    
  
      // local prefix sum for tid
      z[tid] = x[hi-1];
  #pragma omp barrier
  
  
      // global prefix sum over z
      for(j=1; j&lt;nthr; j=2*j) {
        if (tid &gt;= j)
          z[tid] = z[tid] + z[tid-j];
  #pragma omp barrier
      } 
  
      // Update local prefix sum x
      prev_sum = z[tid] - x[hi-1];
      for(i=lo;i&lt;hi;i++)
        x[i] += prev_sum;
    }
    free(z);
  }
  
  void initlist(int *iplist, unsigned long size)
  {
    int i;
    for ( i = 0; i &lt; size; i++)
      iplist[i] = i+1;
    //     iplist[i] = rand() % 13;
  }
  
  void printlist(int *list, unsigned long size)
  {
    int i;
    for(i = 0; i &lt; size; i++) {
      printf("%d ", list[i]);
    }
    printf("\n");
  }
  
  void computeseqprefixsum(int *iplist, int *seqprefixsum, unsigned long size)
  {
    int i;
    seqprefixsum[0] = iplist[0];
    for(i = 1; i &lt; size; i++) {
      seqprefixsum[i] = seqprefixsum[i-1] + iplist[i];
    }
  }
  
  void checkresults(int *seqsum, int *parallelsum, unsigned long size)
  {
    int i;
    for(i = 0; i &lt; size; i++)
    {
      if(seqsum[i] != parallelsum[i]) {
        printf("seqsum[%d] = %d != parallelsum[%d] = %d\n", i, seqsum[i], i,
            parallelsum[i]);
        exit(1);
      }
    }
  }
  
  int main(int argc, char *argv[])
  {
    // seed the rand generator
    srand(time(NULL));
    double seqstart, seqend, parstart, parend, seqtime, partime;
  
    // initialize list
  
    int *iplist, *seqprefixsum, *pprefixsum ;
  
    iplist = (int*) malloc(sizeof(int) * N);
    seqprefixsum = (int*) malloc(sizeof(int) * N);
    pprefixsum = (int*) malloc(sizeof(int) * N);
  
    if(iplist == NULL || seqprefixsum == NULL || pprefixsum == NULL) {
      printf("memory cannot be allocated\n");
      exit(1);
    }
  
    initlist(iplist, N);
  
    seqstart = omp_get_wtime();
  
    computeseqprefixsum(iplist, seqprefixsum, N);
  
    seqend = omp_get_wtime();
  
    seqtime = seqend - seqstart;
  
    omp_set_num_threads(NUM_THREADS);
  
    parstart = omp_get_wtime();
  
    computeparallelprefix(iplist, pprefixsum, N);
  
    parend= omp_get_wtime();
  
    partime = parend - parstart;
  
    checkresults(seqprefixsum, pprefixsum, N);
  
    printf("Seq Time : %f, Par Time : %f, Speedup : %f\n", seqtime, partime,
        seqtime/partime);
  
    free(iplist); free(seqprefixsum); free(pprefixsum);
  
    return 0;
  }
  </code></pre>