//
// Program to solve 1D heat-distribution problem (1D version of
//   problem presented in textbook, section 6.3.2)
//
// Command-line arguments:  
//   number of threads
//   number of interior points (must be a multiple of number
//     of processes)
//   value for left end
//   value for right end
//   maximum iterations (optional)
//   convergence threshold (optional)
//
// Output:
//   status messages to standard error
//   values of all points on convergence to standard output
//
// Version 2:  Starts up P threads, each of which performs
//   the main loop, using barriers to synchronize.
//
#include <iostream>
#include <iomanip>		// has setprecision()
#include <cstdlib>		// has exit(), etc.
#include <algorithm>		// has fill()
#include <cmath>		// has fabs()
#include "timer.h"		// has timer()
#include "pthreads-threadmgr.h"	// has threadManager class
#include "pthreads-barrier.h"	// has barrier class

#define DEFAULT_THRESHOLD      0.01     // convergence threshold
#define DEFAULT_MAX_ITERATIONS 1000

// ---- Class for doing calculations --------------------------------

// This will be an "enclosing class" to hold variables to be shared
//   among threads.  We'll build one object of this class, and pass
//   a pointer to it to all the threads.  This seems tidier than
//   either of the previous approaches (global variables, or static
//   class variables).  This class will also provide methods that
//   do the actual desired work, including starting up threads and
//   waiting for them to complete.

class computeObj {
public:
  // Constructor.
  computeObj(const int nThreads, const int nPoints, 
	     const int maxIter, const double threshold,
	     const double leftEnd, const double rightEnd);
  // Destructor.
  ~computeObj(void);

  // Do calculations.
  void calculate(void);

  // Do output.
  void output(void);

private:
  // These variables hold data to be shared among all threads.
  int nThreads;
  int nPoints;
  int maxIter;			// maximum iterations
  double threshold;		// convergence threshold
  double leftEnd;		// left boundary point
  double rightEnd;		// right boundary point
  int iter;			// number of iterations
  bool converged;		// global convergence?
  double * oldval;		// array of old values
  double * newval;		// array of new values
  bool * convergedLocal;	// per-thread convergences
  barrier * barPtr;

  // ---- Class for thread -------------------------------------------
  class computeThreadObj {
  public:
    typedef computeThreadObj * ptr;
    computeThreadObj(const int myID,
		     const int firstIndex, const int lastIndex,
		     computeObj * context);
    void run(void);
  private:
    int myID;		
    int firstIndex;		// index of first element we "own"
    int lastIndex;		// index of last element we "own"
    computeObj * context;	// pointer to "enclosing object"
  };
};

// ---- Main program -------------------------------------------------

int main(int argc, char* argv[]) {

  // Process command-line arguments.

  if (argc < 5) {
    cerr << "Usage:  " << argv[0] 
	 << " nThreads nPoints leftTemp rightTemp"
	 << " [maxIter] [threshold]\n";
    exit(EXIT_FAILURE);
  }

  int nThreads = atoi(argv[1]);
  int nPoints = atoi(argv[2]);
  if ((nPoints % nThreads) != 0) {
    cerr << "nPoints must be a multiple of nThreads\n";
    exit(EXIT_FAILURE);
  }

  int maxIter = DEFAULT_MAX_ITERATIONS;
  if (argc > 5) 
    maxIter = atoi(argv[5]);

  double threshold = DEFAULT_THRESHOLD;
  if (argc > 6)
    threshold = atof(argv[6]);

  cerr << "Computing for " << nPoints << " points using "
       << nThreads << " threads\n";

  // Build object, using values for boundary points from command line.

  double startInitTime = timer();
  computeObj obj(nThreads, nPoints, maxIter, threshold, 
		 atof(argv[3]), atof(argv[4]));

  // Do calculations.
  double startTime = timer();
  obj.calculate();
  double endTime = timer();

  // Print results.
  obj.output();
  double endOutputTime = timer();
  cerr << "Time for computation (seconds) = " << setprecision(4) 
       << endTime - startTime << endl;
  cerr << "Total time (seconds) = " << setprecision(4)
       << endOutputTime - startInitTime << endl;

  // Exit.
  return EXIT_SUCCESS;
}

// ---- Functions for computeObj class -------------------------------

// Constructor.
computeObj::computeObj(const int nThreads, const int nPoints, 
		       const int maxIter, const double threshold,
		       const double leftEnd, const double rightEnd) {
  this->nThreads = nThreads;
  this->nPoints = nPoints;
  this->maxIter = maxIter;
  this->threshold = threshold;
  this->leftEnd = leftEnd;
  this->rightEnd = rightEnd;
  oldval = new double[nPoints + 2];
  newval = new double[nPoints + 2];
  convergedLocal = new bool[nThreads];
  barPtr = new barrier(nThreads);
}

// Destructor.
computeObj::~computeObj(void) {
  delete [] oldval;
  delete [] newval;
  delete [] convergedLocal;
  delete barPtr;
}

// Do calculations.
void computeObj::calculate(void) {

  // Initialize arrays:
  //   Interior points are initially = 0.
  fill(oldval+1, oldval+nPoints+1, 0.0);
  //   Endpoints' values are from constructor arguments.
  oldval[0] = leftEnd;
  oldval[nPoints+1] = rightEnd;
  //   Endpoints' values don't change, so set new values here too.
  newval[0] = oldval[0];
  newval[nPoints+1] = oldval[nPoints+1];

  converged = false;		// will become true if, at all points,
				//   abs(oldval - newval) <= threshold.
				// convergedLocal[i] is true iff
				//   abs(oldval - newval) <= threshold
				//   for all points owned by thread i

  iter = 0;			// number of iterations so far

  computeThreadObj::ptr * threadObjs = 
    new computeThreadObj::ptr[nThreads]; // parameters for threads

  // Create computeThreadObj objects, one for each thread, to hold
  //   parameters to pass to threads.  
  int nPointsPerThread = nPoints / nThreads;
  for (int i = 0; i < nThreads; ++i)
    threadObjs[i] = new computeThreadObj(i, 
					 i*nPointsPerThread + 1,
					 (i+1)*nPointsPerThread,
					 this);

  // Create and start nThreads new threads, and wait for them to 
  //   complete.
  threadManager<computeThreadObj> tMgr(nThreads, threadObjs);
  tMgr.join();

  // Clean up.
  for (int i = 0; i < nThreads; ++i) 
    delete threadObjs[i];
  delete [] threadObjs;
}

// Do output.
void computeObj::output(void) {
  if (converged)
    cerr << "Converged in " << iter << " iterations";
  else
    cerr << "Failed to converge in " << iter << " iterations";
  cerr << " (convergence threshold " << threshold << ")\n";
  for (int i = 1; i <= nPoints; ++i) 
    cout << "value for point " << i << " = " << oldval[i] << endl;
}

// ---- Functions for computeObj::computeThreadObj class -------------

// Each thread "owns" one segment of arrays oldval and newval,
//   and is responsible for updating their values.

// Constructor.
computeObj::computeThreadObj::computeThreadObj(const int myID,
					       const int firstIndex, 
					       const int lastIndex,
					       computeObj * context) {
  this->myID = myID;
  this->firstIndex = firstIndex;
  this->lastIndex = lastIndex;
  this->context = context;
}

// Member function containing code each thread should execute.
void computeObj::computeThreadObj::run(void) {

  // Performs main loop:
  //
  //   Each time through, we:
  //     Compute new values for points we "own", using values
  //       computed last time.
  //     Check for local convergence.
  //     Synchronize across all threads.
  //     Combine local convergence results into global convergence.
  //     "Copy" old values to new values (by switching pointers).
  //     Synchronize across all threads.
  //
  //   The loop continues until convergence or the maximum number of
  //     iterations is reached.  At the end, one thread updates
  //     variables needed by the calling program.


  // Initialize more local variables (some are just copies of shared
  //   read-only variables, copied for convenience and possible faster
  //   access).
  double * oldVals = context->oldval;
  double * newVals = context->newval;
  int nThreads = context->nThreads;
  int maxIter = context->maxIter;
  double threshold = context->threshold;
  int iterations = 0;

  // Use shared (global) convergence in loop test (initially false).

  while ((iterations < maxIter) && !(context->converged)) {

    // Compute new values for points (from firstIndex to lastIndex).
    for (int i = firstIndex; i <= lastIndex; ++i)
      newVals[i] = 0.5*(oldVals[i-1] + oldVals[i+1]);

    // Check for convergence in this section.
    bool convergedTemp = true;
    for (int i = firstIndex; i <= lastIndex && convergedTemp; ++i)
      if (fabs(oldVals[i] - newVals[i]) > threshold)
	convergedTemp = false;
    // Save result in this thread's element of the shared array.
    context->convergedLocal[myID] = convergedTemp;

    // Synchronize across all threads.
    context->barPtr->wait();

    // Have one thread check for global convergence.
    if (myID == 0) {
      bool convergedTemp = true;
      for (int i = 0; i < nThreads && convergedTemp; ++i)
	if (!(context->convergedLocal[i]))
	  convergedTemp = false;
      // Save result in shared variable.
      context->converged = convergedTemp;
    }

    // "Copy" old values to new values (by switching pointers).
    double * temp = oldVals;
    oldVals = newVals;
    newVals = temp;

    // Increment count.
    ++iterations;

    // Synchronize across all threads.
    context->barPtr->wait();

  } // end of main loop

  // Update shared variables (in "enclosing object", from which they
  //   will later be printed).
  if (myID == 0) {
    context->iter = iterations;
    context->oldval = oldVals;
    context->newval = newVals;
  }
}