/* --------------------------------------------------------------!
! The current code has the purpose:                              !
! (1) use threadprivate variable in the different program units. ! 
!----------------------------------------------------------------!
! Xiaoxu Guan at LSU HPC.                                        !
! March 07, 2016                                                 !
!-------------------------------------------------------------- */

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define nsize 410


/* global variables */
int id, istart,iend;
double array_a[nsize-1],array_b[nsize-1],array_c[nsize-1];
#pragma omp threadprivate (id,istart,iend) 
/* end of global variables */


int main () {

int i,nothread,nchunk;
double diff;
double array_master[nsize-1];


   for(i = 0; i < nsize; i++ ) 
   {
      array_a[i] = (double) i;
      array_b[i] = (double) (i+i);
   }


#pragma omp parallel
    {
      nothread = omp_get_num_threads();
      id = omp_get_thread_num();
      nchunk = (double)(nsize) / (double)(nothread);
      istart = id*nchunk;
      iend = istart + nchunk - 1;
      if(id == nothread-1) {
         istart = id*nchunk;
         iend = nsize-1;
         nchunk = iend - istart + 1;
         printf("Hello from %5d %5d %5d\n", id, istart, iend);
        }
    }


       printf("Hello from master. id = %d\n", id);


#pragma omp parallel 
    {
      scalar();
    }

/* check the results:  */
      for( i = 0; i < nsize; i++) {
         array_master[i] = array_a[i] + array_b[i];
         diff = ( array_master[i] - array_c[i] ) / array_master[i];
         if( fabs(diff) >= 1.e-10 ) {
             printf("Error seems too large!\n");
             printf("i = %d diff = %f %d %d\n", i, diff, array_master[i], array_c[i]);
            } 
      }

}   /* end of main */




/*--------------------------------------------------!
! this routine needs to be called by all threads.   !
!--------------------------------------------------*/

int scalar() {

int i;

/* each thread only works on its own chunk of arrays.  */
   printf("from scalar routine. id = %5d  istart = %5d  iend = %5d\n", id, istart, iend);

   for (i = istart; i <= iend; i++) {
        array_c[i] = array_a[i] + array_b[i];
       }

}


