/****************************************************************************************************
openmp_vnorm.c calculates a vector norm using OpenMP. The program generates a random vector of length specified by the user and calculates the norm. Parallelism is achieved by having each thread calculate the sum of the squares of a section of the vector; the OpenMP Reduction() clause is used to consolidate those values into a single variable (the square of the vector norm).

The program takes two arguments: 
(1) The first is the number of threads to be used. If no arguments are provided, the program assumes 4 threads.
(2) The second is the length of the vector. If only the first argument is provided, the program uses 100.

The program prints the vector norm and writes the vector to the file 'vec.txt'.

Author: Justin Krometis
Modified: February 20, 2012
****************************************************************************************************/

#include <stdio.h>
#include <string.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>

//Functions to be used ( see definitions below main() )
void fprint_vec(char flnm[], double *vec, int dim);
void rand_vec(double *vec, int dim);

//main() does the work
int main(int argc, char *argv[] ) {
  double  *vec,     //vector
          n_sqr=0;  //square of the vector norm
  
  int nthr,   //number of threads
      dim,    //length of the vector
      i;      //iterator

  //the first argument is the number of threads. if none is provided, assume 4.
  if(argc>1)
    nthr=atoi(argv[1]);  
  else {
    //directive in case a user doesn't know the parameters
    printf("\nProgram parameters are: [number of threads] [vector dimension].");
    printf("\nNo number of threads provided. Assuming 4.");
    nthr=4;
  }

  //the second argument is the length of the vector. if none is provided, use 100.
  if(argc>2)
    dim=atoi(argv[2]);
  else {
    printf("\nNo vector dimension provided. Using 100.");
    dim=100;
  }

  //allocate memory for the three matrices
  vec = (double *)malloc(dim*sizeof(double));
  if(vec==NULL) printf("\nUnable to allocate memory for the vector.\n");
  
  //populate the vector with random numbers
  printf("\nGenerating a random vector of length %d...",dim);
  rand_vec(vec, dim);
  printf("Done.\nCalculating the norm...");

  #pragma omp parallel shared(vec, dim) private(i) //create threads
  {
    //Split up the for loop
    //Use the reduction() clause to have OpenMP sum up all of the private copies of n_sqr
    #pragma omp for reduction(+:n_sqr) 
    for(i=0; i<dim; i++) //iterate through the rows of the result
    {
      printf("\nRow %d being handled by thread ID %d.",i,omp_get_thread_num()); //uncomment this line to see which thread handles each row
      n_sqr += vec[i]*vec[i];
    }
  }
  
  printf("Done.\nVector norm is %.5f.\n\n",sqrt(n_sqr));

  fprint_vec("vec.txt",vec,dim);

  free(vec); //free the memory allocated for the vector

  return 0;
} // end main()

//fprint_vec() prints vector 'vec' of length 'dim' to the file 'flnm'
void fprint_vec(char flnm[], double *vec, int dim)
{
  int i;
  FILE *fl;

  fl = fopen(flnm,"w"); //open the file

  for(i=0; i<dim; i++)
    fprintf(fl,"%1.4f ",*( vec+i ));

  fclose(fl); //close the file
} //end fprint_vec()

//rand_vec() generates a random vector of length 'dim'
void rand_vec(double *vec, int dim)
{
  int i; //iterators
 
  srand( (unsigned int)time(NULL) );
  for(i=0; i<dim; i++) //iterate through the vector
    *(vec+i) = (10.0*rand())/RAND_MAX;
} //end rand_vec()