UNB/ CS/ David Bremner/ teaching/ cs3383/ examples/ mvdc.cc
#include <random>
#include <vector>
#include <iostream>
#include <cilk/cilk.h>
#include <time.h>
#include <unistd.h>
#include <climits>
#include "mat.hpp"

using namespace std;

void
RowMult(const mat &A, const vec &x, vec &y, int i) {
  y[i]=0;
  for (int j=0; j < A[i].size(); j++){
    y[i] += A[i][j]*x[j];
  }
}

void MVDC(const mat &A, const vec &x, vec &y, int f, int t) {
  if (f == t) {
    RowMult(A, x, y, f);
  } else {
    int m = (f+t)/2;
    cilk_spawn MVDC(A,x,y,f,m);
    MVDC(A,x,y,m+1,t);
    cilk_sync;
  }
}

void usage (int argc, char** argv) {
  cerr << "usage: " << argv[0] << " [-d] [-s seed] <rows> <cols> " << endl;
  exit(1);
}

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

  int rows,cols;

  int opt;
  int debug = 0, seed_found = 0;
  long seed;
  
  while ((opt = getopt(argc, argv, "ds:")) != -1) {
    switch (opt) {
    case 'd':
      debug = true;
      break;
    case 's':
      seed = atol (optarg);
      seed_found = 1;
      break;
    default:
      usage(argc, argv);
    }
  }
  
  if (optind + 1 >= argc)
    usage(argc, argv);
    
  rows = atol(argv[optind]);
  cols = atol(argv[optind+1]);
  
  if (! seed_found) {
    random_device rd;
    seed = rd();
  }

  clock_t start = clock();
  
  mat A=random_mat(seed, rows, cols);

  vec x,y;

  if (debug)
    cout << "x = " << endl;

  mt19937 gen(seed+rows+cols);
  for(int j=0; j < cols; j++) {
    scalar val = dist(gen());
    x.push_back(val);
    if (debug) {
      cout << " " << val;
    }
  }
  if (debug)
    cout << endl;

  y.resize(rows);

  clock_t init = clock();

  MVDC(A,x,y,0,A.size() - 1);
  
  clock_t end = clock();
  
  cout << "init = " << (double)(init - start)/CLOCKS_PER_SEC << endl << "mult = " << (double)(end - init)/CLOCKS_PER_SEC  << endl;

  if (debug) {
    cout << "y =" << endl;
    for(int j=0; j < rows; j++) {
      cout << " "  << y[j];
    }
    cout << endl;
  }
  
}