///////////////////////////////////////////////////////////////////////
//ompEx.cpp
#include <omp.h>
#include <iostream>

int main () {
  int nP = omp_get_num_procs();
  omp_set_num_threads(nP);
  int myNumber;
#pragma omp parallel private(myNumber)
  {
    myNumber = omp_get_thread_num();
#pragma omp critical
    {
    std::cout << "Hi ya'll from thread " << 
      myNumber << "!" << std::endl;
    }
  }
  return 0;
}
///////////////////////////////////////////////////////////////////////
//ompOptDriver.cpp
#include <iostream>
#include <cmath>
#include <omp.h>
#include "functor.h"
#include "optim.h"

double pi = 2.*acos(0.);
 
using std::cin; using std::cout; using std::endl;

double f(double theta){
  return theta*sin(4*pi*theta*theta);
}

int main(int argc, char** argv){
  Functor func(&f);
  
  int nEvals;
  cin >> nEvals;

  int nP = omp_get_num_procs();
  omp_set_num_threads(nP);

  double myMin; 
  double* pMin = new double[nP];
  double myLow, myUp;
  int myRank;
  Optim Opt(nEvals);
 
#pragma omp parallel private(myRank, myMin, myLow, myUp) 
#pragma omp firstprivate(Opt, func)
  {
    myRank = omp_get_thread_num();
    myLow = (double)myRank/(double)(nP);
    myUp = (double)(myRank + 1)/(double)(nP);
    myMin = Opt.golden(func, myLow, myUp);
    pMin[myRank] = myMin;
  }

  double min = pMin[0];
  cout << min << endl;
  for(int i = 1; i < nP; i++){
    if(func(pMin[i]) < func(min)) min = pMin[i];
    cout << pMin[i] << endl;
  }

  cout << "The minimizer is " << min << endl; 
  return 0;
}

#ompOpt.mk
ompOpt : ompOptDriver.o optim.o ranGen.o
	g++ -fopenmp ompOptDriver.o optim.o ranGen.o -o ompOpt  

optim.o : optim.cpp optim.h  ranGen.h
	g++ -c -fopenmp optim.cpp

ranGen.o : ranGen.cpp ranGen.h
	g++ -c -fopenmp ranGen.cpp

ompOptDriver.o : ompOptDriver.cpp optim.h
	g++ -c -fopenmp ompOptDriver.cpp
///////////////////////////////////////////////////////////////////////
//mpiEx1.cpp                                                               
#include <iostream>
#include <mpi.h>

int main(){

  MPI::Init();

  int myRank = MPI::COMM_WORLD.Get_rank();
  int nP = MPI::COMM_WORLD.Get_size();

  std::cout << "Hello ya'll from processor " << myRank;
  std::cout << " of " << nP << std::endl;

  MPI::Finalize();
  return 0;
}
///////////////////////////////////////////////////////////////////////
//mpiEx2.cpp
#include <iostream>
#include "mpi.h"

int main(){

  MPI::Init();
  int myRank = MPI::COMM_WORLD.Get_rank();
  int nP = MPI::COMM_WORLD.Get_size();

  char msg[20];
  if(myRank == 0){
    strcpy(msg, "Hello ya'll");
    for(int i = 1; i < nP; i++)
      MPI::COMM_WORLD.Send(msg, 20, MPI::CHAR, i, 0);
  }
  else{
    MPI::COMM_WORLD.Recv(msg, 20, MPI::CHAR, 0, 0);
    std::cout << msg << " from processor " << myRank << std::endl;
  }
  
  MPI::Finalize();
  return 0;
}
///////////////////////////////////////////////////////////////////////
//summary.h
#ifndef SUMMARY_H
#define SUMMARY_H
class Summary{
  int n;//sample size and number of variables
  double* pX;//pointer to data array
  
public:
  Summary(){pX = 0;}
  Summary(int N, double* pData){
    n = N;
    pX = pData;
  }	
  ~Summary(){delete[] pX;};//destructor
  void west(double* stat);
};
#endif

//summary.cpp
#include "summary.h"
#include <cmath>

void Summary::west(double* stat){

    double xbar = pX[0];
    double RSS = 0;
    for(int i=0; i< n; i++){
      RSS += ((double)(i)/(double)(i+1))*(pX[i] - xbar)*(pX[i] - xbar);
      xbar += (1./(double)(i+1))*(pX[i] - xbar);
    }
    stat[0] = xbar;
    stat[1] = RSS/(double)(n - 1);
}

//anovaMPI.cpp
#include <iostream>
#include <fstream>
#include <cstdlib>
#include "mpi.h"
#include "summary.h"

using std::cout; using std::endl;
using std::ifstream;

int main(int argc, char** argv){
  double** pData;
  double* pMyData;

  MPI::Init(argc, argv);
  
  int myRank = MPI::COMM_WORLD.Get_rank();
  int nP = MPI::COMM_WORLD.Get_size();
  int n = atoi(argv[1]);

  if(myRank == 0){
    pData = new double*[nP];
 
    ifstream fIn;
    fIn.open(argv[2]);
    for(int j = 0; j < nP; j++){
      pData[j] = new double[n];
      for(int i = 0; i < n; i++){
	fIn >> pData[j][i];
      }
    }

    fIn.close();
      
    for(int i = 1; i < nP; i++)
      MPI::COMM_WORLD.Send(pData[i], n, MPI::DOUBLE, i, 0);
    Summary mySummary(n, pData[0]);
    double* stat = new double[2];
    mySummary.west(stat);
    cout << "For processor " << myRank 
	 << " the mean and variance are " 
	 << stat[0] << " and " << stat[1] << endl;

    double num = 0;
    double xBar = stat[0];
    double den = ((double)(n - 1))*stat[1];
    
    for(int i = 1; i < nP; i++){
      MPI::COMM_WORLD.Recv(stat, 2, MPI::DOUBLE, i, 1);
      num +=  ((double)(i)/(double)(i + 1))*(stat[0] - xBar)*
	(stat[0] - xBar);
      xBar = ((double)(i)/(double)(i + 1))*xBar + 
	stat[0]/((double)(i + 1));
      den += ((double)(n - 1))*stat[1];
    }
    num *=(double)n/((double)(nP-1));
    den /= (double)(nP*n - nP); 
    cout << "The grand average is " << xBar << endl;
    cout << "The numerator mean square is " << num << " on " 
	 << (nP - 1) << " degrees-of-freedom" << endl;
    cout << "The denominator mean square is " << den << " on " 
	 << (n*nP - nP) << " degrees-of-freedom" << endl;
    if(den != 0)
      cout << "The F-ratio is " << (num/den) << endl;
  }

  else{
    double* mypData = new double[n];
    MPI::COMM_WORLD.Recv(mypData, n, MPI::DOUBLE, 0, 0);
    Summary mySummary(n, mypData);
    double* stat = new double[2];
    mySummary.west(stat);
    cout << "For processor " << myRank 
	 << " the mean and variance are " 
	 << stat[0] << " and " << stat[1] << endl;
    MPI::COMM_WORLD.Send(stat, 2, MPI::DOUBLE, 0, 1);
    }  

  MPI::Finalize();

  return 0;
}
///////////////////////////////////////////////////////////////////////
//ranOpenMP.cpp
#include <omp.h>
#include "RngStream.h"
#include <iostream>
using std::cout; using std::endl; using std::cin;

int main(){
  int nP = 4;
  omp_set_num_threads(nP);                              
  unsigned long seed[6] ={1,2,3,4,5,6};

  RngStream::SetPackageSeed (seed);
  RngStream RngArray[nP];
  for(int i = 0; i < nP; i++){
    cout << "For array element " << i << endl;
    RngArray[i].WriteState();
    cout << "The random number is " 
         << RngArray[i].RandU01() << endl;
  }

  int myRank;
  RngStream Rng;
  RngStream::SetPackageSeed (seed);

#pragma omp parallel private(Rng, myRank)
  {
  myRank = omp_get_thread_num();
#pragma omp critical
    {
      cout << "For thread " <<  omp_get_thread_num() << endl;
      Rng.WriteState();
      cout << "The random number is " << Rng.RandU01() << endl;
    }
  }
  return 0;
}
///////////////////////////////////////////////////////////////////////
//ranMPI.cpp
#include <iostream>
#include <fstream>
#include "mpi.h"
#include "RngStream.h"

using std::cout; using std::endl;

int main(){
  unsigned long seed[6] ={1,2,3,4,5,6};
 
  MPI::Init();
  int myRank = MPI::COMM_WORLD.Get_rank();
  int nP = MPI::COMM_WORLD.Get_size();

  if(myRank == 0){
    unsigned long tempSeed[6];
    double ranValue = 0;

    //start with the 0 processor
    RngStream::SetPackageSeed(seed);
    RngStream Rng;

    double myNumber = Rng.RandU01();
    cout << "For processor 0 the state is" << endl;
    for(int j = 0; j < 6; j++)
      cout << seed[j] << endl;
    cout << endl;

    //create and send the seeds to the other processors    
    for(int i = 1; i < nP; i++){
      RngStream::AdvanceSeed(seed, tempSeed);
      MPI::COMM_WORLD.Send(tempSeed, 6, MPI::UNSIGNED_LONG, i, 0);
      cout << "For processor " << i << " the state is" << endl;
      for(int j = 0; j < 6; j++){
	cout << tempSeed[j] << endl;
      }
      cout << endl;
      for(int j = 0; j < 6; j++)
	seed[j] = tempSeed[j];
    }
     
  //now collect the results
    for(int i = 1; i < nP; i++){
      MPI::COMM_WORLD.Recv(&ranValue, 1, MPI::DOUBLE, i, 0);
      cout << "The random number for processor " << i 
	   << " is "<< ranValue << endl;  
    }
    cout << "The random number for process 0 is " << myNumber
	 << endl;
  }   
  
  else{
    unsigned long mySeed[6];
    //get my seed
    MPI::COMM_WORLD.Recv(mySeed, 6, MPI::UNSIGNED_LONG, 0, 0);

    //create the rngstream object
    RngStream::SetPackageSeed(mySeed);
    RngStream Rng;
    
    //generate my number and send it back to the 0 processor
    double myNumber = Rng.RandU01();
    MPI::COMM_WORLD.Send(&myNumber, 1, MPI::DOUBLE, 0, 0);
  }

  MPI::Finalize();

  return 0;
}
///////////////////////////////////////////////////////////////////////