prefix/hillis_sum.c

   1 /*
   2  * O(nlog n) work algorithm (Hillis-Steele)
   3  */
   4 #include "hillis_sum.h"
   5
   6 /*
   7  * Hillis/Steele, prefix sum version
   8  */
   9 unsigned long is_power_of_two(unsigned long var) {
  10         if (var == 1) {
  11                 return 1;
  12         }
  13         while (var > 2) {
  14                 if ((var & 0x1) == 0x1) {
  15                         return 0;
  16                 }
  17                 var >>= 1;
  18         }
  19         return 1;
  20 }
  21
  22 void algorithm (numtype x[], unsigned long size, unsigned int ops[]) {
  23         unsigned long k;
  24         unsigned long i;
  25
  26         if (is_power_of_two(size) != 1) {
  27                 fprintf(stderr, "Hillis-Sum algorithm only works on arrays on size 2^n\n");
  28                 exit(1);
  29         }
  30
  31         for(k=2; k <= size; k <<=1){
  32                 #pragma omp parallel for shared(x, size, ops, k) private(i)
  33                 for(i = (k-1); i < size; i+=k){
  34 #ifdef DEBUG
  35                         printf ("x[%2li] = x[%2li] + x[%2li]; // {i:%li, k:%li}\n", i, i-k, i, i, k);
  36 #endif
  37                         x[i] = x[i-(k/2)] + x[i];
  38                 }
  39         }
  40 }
  41