openmp/prefix/hillis.c

   1 /*
   2  * O(nlog n) work algorithm (Hillis-Steele)
   3  */
   4 #include "hillis.h"
   5
   6 /*
   7  * Hillis/Steele, prefix sum version
   8  */
   9 /*
  10 void algorithm (numtype x[], unsigned long size, unsigned int ops[]) {
  11         numtype tmp[size];
  12         unsigned long k;
  13         unsigned long i;
  14
  15         for(k=2; k <= size; k <<=1){
  16                 #pragma omp parallel for shared(x, size, ops, k) private(i)
  17                 for(i = (k-1); i < size; i+=k){
  18 #ifdef DEBUG
  19                         printf ("x[%2li] = x[%2li] + x[%2li]; // {i:%li, k:%li}\n", i, i-k, i, i, k);
  20 #endif
  21                         x[i] = x[i-(k/2)] + x[i];
  22                 }
  23         }
  24 }
  25 */
  26
  27 /*
  28  * Hillis/Steele, partial prefix sum version
  29  */
  30 void algorithm (numtype x[], unsigned long size, unsigned int ops[]) {
  31         numtype *y;
  32         unsigned long k;
  33         unsigned long i;
  34
  35         y = malloc(size * sizeof(numtype));
  36         y[0] = x[0];
  37
  38         for (k=1; k<size; k<<=1) {
  39                 #pragma omp parallel for shared (x, y, ops, k) private (i)
  40                 for (i=k; i<size; i++) {
  41 #ifdef DEBUG
  42                         printf ("y[%li] = x[%li] + x[%li]; // %li + %li", i, i-k, i, x[i-k], x[i]);
  43 #endif
  44                         y[i] = x[i-k]+x[i];
  45                 }
  46                 memcpy (x, y, size * sizeof(numtype));
  47         }
  48
  49         if (y != NULL) {
  50                 free(y);
  51                 y = NULL;
  52         }else{
  53                 printf ("Has not been freed!!!\n");
  54         }
  55 }