Intel Cilk Plus -- Example #5 (fill the matrix and calculate column prefix sums)

fill.cpp

#include <cstdlib>
#include <cilk/cilk.h>
#include <cilk/reducer.h>
#include <cilk/reducer_opadd.h>
#include <iostream>
#include <assert.h>
#include <time.h>
#include <sys/time.h>

/* ex: set tabstop=8 softtabstop=4 expandtab: */

void prng_power(int &, int &, int, int);
int prng_jump(int, int, int, int, int);
static inline unsigned long long cilk_getticks();

static int * fill_matrix(int M, int N, int seed0, int a, int b, int m, int direction)
{
    int area = N * M;
    int ax = 1, bx = 0;
    int ay = 1, by = 0;

    unsigned long long start_tick, end_tick;
    start_tick = cilk_getticks();

    seed0 = (a * seed0 + b) % m;
    if (direction == 1) {
	ax = a; bx = b;
	ay = a; by = b;
	prng_power(ay, by, m, N);
    } else {
	ax = a; bx = b;
	ay = a; by = b;
	prng_power(ax, bx, m, N);
    }

    // fill the first row
    int * row = (int *) _mm_malloc(N * sizeof(int), 64);
    {
        int seed = seed0;
        for(int j = 0; j < N; j++) {
            row[j] = seed;
            seed = (ax * seed + bx) % m;
        }
    }

    // calculate the mean PRNG value
    cilk::reducer_opadd<long long> sum;
    int mean, remainder;
    cilk_for(int i = 0; i < M; i++) {
        int ay2 = ay;
        int by2 = by;
        prng_power(ay2, by2, m, i);
	sum += __sec_reduce_add((ay2 * row[0:N] + by2) % m);
    }
    mean = (int) (sum.get_value() / (long long) area);
    remainder = (int) (sum.get_value() % (long long) area);
    mean += (remainder * 2 > (signed) area) ? (1) : (0);

    // fill matrix with PRNG values, corrected by the mean value
    int * A = (int *) _mm_malloc(N * M * sizeof(int), 64);
    assert(A != NULL);
    cilk_for(int i = 0; i < M; i++)
    {
        int ay2 = ay;
        int by2 = by;
        prng_power(ay2, by2, m, i);
        int *A1 = (int *) &(A[i * N]);
	A1[0:N] = ((ay2 * row[0:N] + by2) % m) - mean;
    }

    // calculate prefix sums
    for(int i = 1; i < M; i++) {
        int *A0 = (int *) &(A[(i - 1) * N]);
        int *A1 = (int *) &(A[i * N]);
        A1[0:N] += A0[0:N];
    }

    end_tick = cilk_getticks();

    std::cout << "Elapsed time: " << ((int)(end_tick - start_tick)) / 1000.f << "ms" << std::endl;
    std::cout << "PRNG(seed0,a,b,m) sum value: " << sum.get_value() << ", mean value: " << mean << std::endl;

    _mm_free(row);
    return A;
}

int main(int argc, char* argv[])
{
    // int * A = fill_matrix(4, 3, 10, 3, 4, 17, 1);
    // int * A = fill_matrix(1059, 1065, 62, 97, 4, 364, 1);
    int * A = fill_matrix(4096, 65536, 62, 97, 4, 364, 1);
    _mm_free(A);
}

static inline unsigned long long cilk_getticks()
{
     struct timeval t;
     gettimeofday(&t, 0);
     return t.tv_sec * 1000000ULL + t.tv_usec;
}

////// Linear Recurrence Reducer, x[i+1] => (a * x[i] + b) % m
class Recurrence_Reducer {
public:
    static int m;

    // Per-strand view of the data
    struct View {
	friend class Recurrence_Reducer;

    public:
	// Identity value for reducer: a = 1, b = 0
	View() : a(1), b(0) { }

    private:
	int a; 
	int b;
    };

public:
    // View of reducer data
    struct Monoid: cilk::monoid_base<View> {
	static void reduce (View *left, View *right) {
            left->a = (right->a * left->a) % m;
            left->b = (right->b + left->b * right->a) % m;
	}
    };

private:
    // Hyperobject to serve up views
    cilk::reducer<Monoid> imp_;

public:
    Recurrence_Reducer() : imp_() { }

    // Update operations
    inline Recurrence_Reducer& call_next(int a, int b) {
	View &v = imp_.view();
	v.a = (a * v.a) % m;
	v.b = (b + a * v.b) % m; 
	return *this; 
    }
    int get_a() { return imp_.view().a; }
    int get_b() { return imp_.view().b; }
};

int Recurrence_Reducer::m = 0;

void prng_power(int &a, int &b, int m, int n) {
    Recurrence_Reducer rr;
    Recurrence_Reducer s;
    Recurrence_Reducer::m = m;

    // make logarithmic scale
    int p = 0;
    int aa[32], bb[32];

    s.call_next(a, b);
    while (n > 0) {
        int a1 = s.get_a();
        int b1 = s.get_b();
        if (n % 2) {
            aa[p] = a1;
            bb[p] = b1;
        } else {
            aa[p] = 1;
            bb[p] = 0;
        }
        s.call_next(a1, b1);
        n = n >> 1;
        p++;
    }

    // use this scale and calculate n-th power of the linear recurrence
    cilk_for(int i = 0; i < p; ++i) {
        rr.call_next(aa[i], bb[i]);
    }
    a = rr.get_a();
    b = rr.get_b();
}

int prng_jump(int seed0, int a, int b, int m, int n) {
    prng_power(a, b, m, n);
    return (a * seed0 + b) % m;
}
////// END Linear Recurrence Reducer