#include <omp.h>
#include <stdlib.h>
#include <stdio.h>

#define SIZE 100000

int foo(int i)
{
  int result;
  result = (i * i) + 2 * i + 5;
  return result;
}

void copy(int * dst, int * src, int size)
{
  int i;
  for (i = 0; i < size; i++) {
    dst[i] = src[i];
  }
}

int compare(int * a1, int * a2, int size)
{
  int i;
  for (i = 0; i < size; i++) {
    if (a1[i] != a2[i])
      return i;
  }
  return size;
}

int main(int argc, char ** argv)
{
  int a[SIZE];
  int a_seq[SIZE];
  int a_orig[SIZE];
  int b[SIZE];
  int i;
  int n = SIZE;
  int rc;

  for (i=0; i<SIZE; i++) {
    a_orig[i] = rand();
    b[i] = rand();
  }

  // (a)

  for (i = 0; i < n; i++) {
    a_seq[i] = foo(i);
  }

#pragma omp parallel for
  for (i = 0; i < n; i++) {
    a[i] = foo(i);
  }
  if ((rc = compare(a, a_seq, n)) != n) {
    printf("(a) Arrays did not compare at index %d\n", rc);
  } else {
    printf("(a) Arrays compared\n");
  }

  // (b)

  for (i=0; i < n; i++) {
    a_seq[i] = foo(i);
    if (a_seq[i] < b[i]) a_seq[i] = b[i];
  }

#pragma omp parallel for
  for (i=0; i < n; i++) {
    a[i] = foo(i);
    if (a[i] < b[i]) a[i] = b[i];
  }
  if ((rc = compare(a, a_seq, n)) != n) {
    printf("(b) Arrays did not compare at index %d\n", rc);
  } else {
    printf("(b) Arrays compared\n");
  }


  // (c)

  copy(a_seq, a_orig, n);

  int dotp_seq = 0;
  for (i = 0; i < n; i++) {
    dotp_seq += a_seq[i] * b[i];
  }

  // Incorrect solution

  copy(a, a_orig, n);

  int dotp = 0;
#pragma omp parallel for
  for (i=0; i < n; i++) {
    dotp += a[i] * b[i];
  }
  if (dotp != dotp_seq) {
    printf("(c.1) Dot products did not compare\n");
  } else {
    printf("(c.1) Dot products compared\n");
  }

  // Correct solution #1 (but slow)

  int product;

  dotp = 0;
#pragma omp parallel for private(product)
  for (i=0; i < n; i++) {
    product = a[i] * b[i];
#pragma omp critical
    dotp += product;
  }
  if (dotp != dotp_seq) {
    printf("(c.2) Dot products did not compare\n");
  } else {
    printf("(c.2) Dot products compared\n");
  }

  // Correct solution #2

  dotp = 0;
#pragma omp parallel for private(product) reduction(+:dotp)
  for (i=0; i < n; i++) {
    product = a[i] * b[i];
    dotp += product;
  }
  if (dotp != dotp_seq) {
    printf("(c.3) Dot products did not compare\n");
  } else {
    printf("(c.3) Dot products compared\n");
  }

  

  // (d)
  
  copy(a_seq, a_orig, n);

  int k=SIZE/2;

  for (i=k; i< 2*k; i++) {
    a_seq[i] = a_seq[i] + a_seq[i-k];
  }

  copy(a, a_orig, n);

#pragma omp parallel for
  for (i=k; i< 2*k; i++) {
    a[i] = a[i] + a[i-k];
  }

  if ((rc = compare(a, a_seq, n)) != n) {
    printf("(d) Arrays did not compare at index %d\n", rc);
  } else {
    printf("(d) Arrays compared\n");
  }

  // (e)

  int flag;

  copy(a_seq, a_orig, n);
  flag = 0;
  for (i=0; (i < n) && !flag; i++) {
    a_seq[i] = 2.3 * i;
    if (a_seq[i] < b[i]) flag = 1;
  }

  // First note that above formulation does not allow the compiler
  // to determine the number of iterations, and is disallowed by
  // OpenMP (at compile time).  This is the same requirement that
  // disallows a break or goto within the loop.

  // If, however, we reformulate so that we no-op subsequent iterations
  // of the loop, we can try some of the other solutions suggested
  // on the midterm.

  copy(a, a_orig, n);
  flag = 0;
#pragma omp parallel for shared(flag)
  for (i=0; i < n; i++) {
    if (flag) continue;
    a[i] = 2.3 * i;
    if (a[i] < b[i]) flag = 1;
  }
  if ((rc = compare(a, a_seq, n)) != n) {
    printf("(e.1) Arrays did not compare at index %d\n", rc);
  } else {
    printf("(e.1) Arrays compared\n");
  }

  copy(a, a_orig, n);
  flag = 0;
#pragma omp parallel for private(flag)
  for (i=0; i < n; i++) {
    if (flag) continue;
    a[i] = 2.3 * i;
    if (a[i] < b[i]) flag = 1;
  }
  if ((rc = compare(a, a_seq, n)) != n) {
    printf("(e.2) Arrays did not compare at index %d\n", rc);
  } else {
    printf("(e.2) Arrays compared\n");
  }
}