/* File:     omp_dot.c 
 * Purpose:  compute a dot product using OpenMP
 *
 * Input:    n: order of vectors
 *           x, y:  the vectors
 *
 * Output:   the dot product of x and y.
 *
 * Compile:  gcc -g -Wall -fopenmp -o omp_dot omp_dot.c
 * Run:      ./omp_dot <number of threads>
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

void Read_vector(char* prompt, float v[], int n);
float Par_dot(float x[], float y[], int n);
    
int main(int argc, char* argv[]) {
    float* x;
    float* y;
    int    n, thread_count;
    float  dot;

    if (argc != 2) {
       fprintf(stderr, "usage: %s <thread count>\n", argv[0]);
       exit(0);
    }
    thread_count = strtol(argv[1], NULL, 10);

    printf("Enter the order of the vectors\n");
    scanf("%d", &n);
    x = malloc(n*sizeof(float));
    y = malloc(n*sizeof(float));
    Read_vector("the first vector", x, n);
    Read_vector("the second vector", y, n);

    dot = Par_dot(x, y, n);
    printf("The dot product is %f\n", dot);

    free(x);
    free(y);
    return 0;
}  /* main */
   

/*---------------------------------------------------------------
 * Function:  Read_vector
 * Purpose:   Read a vector from stdin
 * In arg:    n:  number of components
 * Out arg:   v:  the vector
 */
void Read_vector(
         char*  prompt  /* in  */,
         float  v[]     /* out */,
         int    n       /* in  */) {
    int i;

    printf("Enter %s\n", prompt);
    for (i = 0; i < n; i++)
        scanf("%f", &v[i]);
}  /* Read_vector */


/*---------------------------------------------------------------
 * Function:  Par_dot
 * Purpose:   Find the dot product of two vectors
 * In args:   x, y:  the two vectors
 *            n:  the number of components in the two vectors
 *            thread_count:  the number of threads
 * Ret val:   the dot product of x and y
 */
float Par_dot(
          float  x[]           /* in */, 
          float  y[]           /* in */, 
          int    n             /* in */) {

    int    i; 
    float  sum = 0.0;
    int    thread_count = omp_get_num_threads();

#   pragma omp parallel for num_threads(thread_count) \
       default(none) shared(n,x,y) private(i) reduction(+:sum)
    for (i = 0; i < n; i++)
        sum += x[i]*y[i];
    return sum;
} /* Par_dot */