// This may look like C code, but it is really -*- C++ -*-

/*
 *  Copyright (c) University of Aizu
 */

#include "Eigen.h"

#ifdef __GNUG__
#pragma implementation
#endif

const int maxdim = 100;

extern "C" { int rand(); }

int randint(int u)
{
  int r = rand();
  if (r < 0) r = -r;
  return r%u;
}

Eigen::Eigen(int dimension)
{
  eigenValues  = new Cvector(n = dimension);
  eigenVectors = new Cmatrix(n, n);
}

Eigen::~Eigen()
{
  delete eigenValues;
  delete eigenVectors;
}

int Eigen::init(int real)
{
  // Create synonymes
  Cvector& lambda = *eigenValues;
  Cmatrix& v      = *eigenVectors;

  // Initial eigenvalues are 0.
  lambda = complex(0);

  // Initial eigenvectors are the canonical basis vectors ...
  v      = complex(1);

  // Solving a simple linear equation should be easy...
  if (n < 2) return 1;

  // if only real eigenvalues expected, don't add spice
  if (real) return 1;

  // ... with a little spice added
  v[n-1][n-2] = v[n-2][n-1] = complex(0,1);

  for (int i = 0; i < n-2; i += 2)
    v[i][i+1] = v[i+1][i] = complex(0,1);

  return 1;
}

int Eigen::print()
{
  for (int i = 0; i < n; i++)
    cout << "> l[" << i << "] = " << (*eigenValues)[i] << endl
	 << "  v[" << i << "] = " << (*eigenVectors)[i];
  return 1;
}

int Eigen::check(const Cmatrix& A, double epsilon, int verbose)
{
  // Create synonymes
  Cvector& lambda = *eigenValues;
  Cmatrix& v      = *eigenVectors;

  double delta;
  int result = 1;

  for (int i = 0; i < n; i++) {
    delta = norm(A * v[i] - lambda[i] * v[i]);
    if (verbose)
      cout << "delta[" << i << "] = " << delta << endl;
    if (delta > n * epsilon)
      result = 0;
  }
  return result;
}

int Eigen::solve(const Cmatrix& A,     // input: matrix
		 double epsilon,       // input: required precision
		 int iterations,       // input: max. number of iterations
		 int randomize,        // input: random scheduling
		 int verbose)          // optional: print intermediary values
{
  // general-purpose indices
  int i,j;

  // Consistancy-checking
  if (n != A.cols()) return 0;
  if (n != A.rows()) return 0;

  // Define processors
  class Processor {
   private:
    int nr;         // number of the processor = dimension assigned
    int done;
    double delta_sum;
   public:

    // Initialise before iterating
    void dimension(int number) { nr = number; done = 0; }
    int startup()              { delta_sum = 0; return !done; }

    // Do one iteration
    void iterate(int j, Cvector& lambda, Cmatrix& v, const Cmatrix &A)
    {
      complex alpha;
      complex delta;

      cout << "+ Processor " << nr << " working on " << j << endl;

      // Step 1
      v[nr] = (A - lambda[j]) * v[nr];

      // Step 2
      alpha = max(v[nr]);

      // Step 3
      v[nr] /= alpha;

      // Step 4
      delta = lambda[nr];
      lambda[nr] = lambda[j] + alpha;
      delta -= lambda[nr];

      // Sum up differences to determine cut-off time
      delta_sum += abs(delta);
    }

    // Clean up after iterating
    void cleanup(double epsilon)
    {
      done = (abs(delta_sum) < epsilon);
      cout << "= Processor " << nr << "delta = " << abs(delta_sum);
      if (done)
	cout << " DONE!";
      cout << endl;
    }
  };
  // Create Processors
  Processor p[maxdim];

  // assign dimensions to each processor
  for (i = 0; i < n; i++)
    p[i].dimension(i);

  // Define schedule
  struct Schedule{
    int processor_nr;
    int iteration_nr;
  };

  // Create Schedule (n processors do n-1 iterations)
  struct Schedule schedule[maxdim*maxdim - maxdim];
  struct Schedule tmp;
  int schedule_size;

  // Do Iterations
  for (; iterations >= 0; iterations--) {
    cout << "@ Iteration " << iterations << endl;
    // create schedule for this round.
    schedule_size = 0;
    // for all processors
    for (i = 0; i < n; i++) {

      // initialize
      if (p[i].startup()) {

	// if work left do be done, insert processor into the schedule
	for (j = 0; j < n; j++)
	  if (i != j) {
	    schedule[schedule_size].processor_nr = i;
	    schedule[schedule_size].iteration_nr = j;
	    schedule_size++;
	  }
      }
    }

    // exit loop if nothing to do.
    if (schedule_size == 0) break;

    // Now randomize schedule
    if (randomize)
      for (i = schedule_size - 1; i > 1; i--) {
	j = randint(i);
	if (j != i) 
	  {tmp = schedule[i]; schedule[i] = schedule[j]; schedule[j] = tmp;}
      }

    // Execute schedule
    for (i = 0; i < schedule_size; i++)
      p[schedule[i].processor_nr].iterate(schedule[i].iteration_nr,
					  *eigenValues,
					  *eigenVectors, A);
    // clean up processors
    for (i = 0; i < n; i++)
      p[i].cleanup(epsilon);
  }

  cout << endl;
  return 1;
}