//
// 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
//
#include <iostream.h>
#include <iomanip>		// has setprecision()
#include <stdlib.h>		// has exit(), etc.
#include <pthread.h>		// has pthread_ routines
#include <algorithm>		// has fill()
#include "threads-timer.h"	// has timer()
#include "threads-threadWrapper.h"	// has threadWrapper class

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

// ---- Class for thread ---------------------------------------------

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

class computeThreadObj {
public:
  typedef computeThreadObj * ptr;

  // member variables are parameters for thread
  int firstIndex;		// index of first element we "own"
  int lastIndex;		// index of last element we "own"
  double threshold;		// convergence threshold
  bool * convergedLocal;	// convergence for this thread's 
				//   segment of the array?
  double * oldVals;		// whole "old values" array
  double * newVals;		// whole "new values" array
  computeThreadObj(const int firstIndex, const int lastIndex,
		   const double threshold, 
		   bool * const convergedLocal,
		   double * const oldVals, double * const newVals)
  {
    this->firstIndex = firstIndex;
    this->lastIndex = lastIndex;
    this->threshold = threshold;
    this->convergedLocal = convergedLocal;
    this->oldVals = oldVals;
    this->newVals = newVals;
  }

  // this member function is the code the thread should execute
  void run(void) {
    // Computes new values for "its" points and determines whether
    //   convergence has occurred.  The convergence variable is
    //   actually this thread's element of an array, so not shared
    //   with other concurrently-executing threads.

    // 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.
    *convergedLocal = true;
    for (int i = firstIndex; i <= lastIndex && *convergedLocal; ++i)
      if (fabs(oldVals[i] - newVals[i]) > threshold)
        *convergedLocal = false;
  }
};

// ---- 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);
  }

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

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

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

  // Set up to time things.

  double startTime = timer();

  // Declare and initialize other variables.

  double * oldval = new double[nPoints + 2];
				// old values
  double * newval = new double[nPoints + 2];
				// new values

  // interior points initially = 0
  fill(oldval+1, oldval+nPoints+1, 0.0);

  // boundary points' values given by command-line arguments
  oldval[0] = atof(argv[3]);
  oldval[nPoints+1] = atof(argv[4]);

  // boundary points do not change, so set new values here
  newval[0] = oldval[0];
  newval[nPoints+1] = oldval[nPoints+1];

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

  unsigned int iter = 0;	// number of iterations so far

  computeThreadObj::ptr * threadObj = 
    new computeThreadObj::ptr[nThreads]; 
				// parameters for threads
  pthread_t * threads = new pthread_t[nThreads];
				// "handles" 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)
    threadObj[i] = new computeThreadObj(i*nPointsPerThread + 1,
					(i+1)*nPointsPerThread,
					threshold, 
					&convergedLocal[i],
					oldval, newval);

  // Main processing loop:
  //   Each time through, we:
  //     Compute new values for all interior points, using values
  //       computed last time.
  //     Check for convergence.
  //     "Copy" new values to old values (by switching pointers).
  //   The loop continues until convergence or the maximum number of
  //     iterations is reached.
  //   Concurrency is obtained by dividing up the work of computing
  //     new values and checking for convergence among nThreads threads.

  for ( ; (iter < maxIterations) && !converged; ++iter) {

    // Start nThreads new threads, each running wrapper function,
    //   with a computeThreadObj object as parameter.  (The wrapper
    //   function just invokes the object's "run()" method.)

    for (int i = 0; i < nThreads; ++i) {
      // update parameters
      threadObj[i]->oldVals = oldval;
      threadObj[i]->newVals = newval;
      if (pthread_create(&threads[i], NULL, 
			 threadWrapper<computeThreadObj::ptr>::run,
			 (void *) threadObj[i]) != 0)
	cerr << "Unable to create thread " << i << endl;
    }

    // Wait for all threads to complete.

    for (int i = 0; i < nThreads; ++i) {
      if (pthread_join(threads[i], NULL) != 0)
	cerr << "Unable to perform join on thread " << i << endl;
    }

    // "Copy" old values to new values by switching pointers.

    double * temp = oldval;
    oldval = newval;
    newval = temp;

    // Check for overall convergence.

    converged = true;
    for (int i = 0; i < nThreads && converged; ++i) 
      converged = converged && convergedLocal[i];

  } // end of main loop
 
  // Print results.

  double endTime = timer();
  cerr << "Elapsed time (seconds) = " << setprecision(4) 
       << endTime - startTime << endl;
  
  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;

  // Free everything allocated with "new" in this function.

  delete [] oldval;
  delete [] newval;
  delete [] convergedLocal;
  for (int i = 0; i < nThreads; ++i)
    delete threadObj[i];
  delete [] threadObj;
  delete [] threads;

  return EXIT_SUCCESS;
}