// tsnran1.cc
// Nick Macks (macks@fnal.gov), Summer 2002

// LinearAlgebra: invert a random 5x5 or a random 6x6 Matrix
// (user can specify if the matrix is to be declared specialized)
// (user can specify any dimension from 1 to 6, not just 5 or 6)

#include "LinearAlgebra/Matrix.h"
#include "ZMtools/ZMtimer.h"
#include "CLHEP/Random/MTwistEngine.h"
#include "CLHEP/Random/RandomEngine.h"

using namespace std;

// Function declarations
void Identification();
void GetDimension(int&);
void Describe(const int&);
void BuildMatrix(Matrix&, const int&);
void DeclareMatrixSpecialized(Matrix&, const int&);
void Show(const Matrix&);
void GetCalcs(int&);
void GetRuns(int&);
void Invert(Matrix&, const int&, const int&, double&, double&);
void Summary(const double&, const double&);


// Package ID
void Identification()
{
  cout << "Using LinearAlgebra\n";
} 

// Ask for the desired dimension
void GetDimension(int& dimension)
{
  cout << "Enter the desired dimension(1-6): ";
  cin >> dimension;  

  while ( (dimension < 1) || (dimension > 6) )
  {
    cout << "Enter an integer between 1 and 6: ";
    cin >> dimension;
  }
}

// Describe how the matrix was defined
void Describe(const int& dimension)
{
  cout << "Matrix was defined 'Matrix ("
       << dimension << ", " << dimension << ")'\n";
}

// Build the matrix
void BuildMatrix(Matrix& matrix, const int& dimension)
{
  // Declare a random number generating engine
  MTwistEngine my_engine(2);

  for (int i=0; i<dimension; ++i)
  {
    for (int j=0; j<dimension; ++j)
    {
      // Get a random number of type double between 0 and 1
      double rand = my_engine.flat();

      // Assign it an integer between 1 and 10
      matrix(i,j) = int(rand*10) + 1;
    }
  }
}

// Ask the user if he wants the matrix declared specialized
void DeclareMatrixSpecialized(Matrix& matrix, const int& dimension)
{
  char special;
  cout << "Declare Matrix as Specialized? (y/n): ";
  cin >> special;

  while ( (special != 'y') && (special != 'Y')
       && (special != 'n') && (special != 'N') )
  {
    cout << "Please enter y or n: ";
    cin >> special;
  }
  
  if ( (special == 'y') || (special == 'Y') )
  {
    if (dimension == 5)
      matrix.declareM55();
    else if (dimension == 6)
      matrix.declareM66();
    else if (dimension == 4)
      matrix.declareM44();
    else if (dimension == 3)
      matrix.declareM33();
    else if (dimension == 2)
      matrix.declareM22();
    else if (dimension == 1)
      cout << "Matrix wasn't declared as specialized since dimension is 1\n";
  }
}

// Show the matrix
void Show(const Matrix& matrix)
{
  cout << "Using matrix =\n" << matrix << "\n";
}

// Get the number of inversions
void GetCalcs(int& N)
{
  cout << "How many inversions? ";
  cin >> N;
}

// Get the number of runs
void GetRuns(int& run)
{
  cout << "How many runs? ";
  cin >> run;
  cout << "\n";
}

// Invert the matrix
void Calculate(Matrix& matrix, const int& N, const int& runs,
               double& minimum, double& maximum)
{
  // Define a timer
  ZMcpuTimer timer;

  // Define a counter
  int counter = 0;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Start the timer
    timer.start();

    // Do the inversions
    for (int i=0; i<N; ++i)
    {
      matrix.inverse();
      ++counter;
    }

    // Stop the timer
    timer.stop();

    // Assign mininum and maximum values
    // on the first run of the loop
    if (r == 0)
    {
      minimum = timer.cpuTime();
      maximum = timer.cpuTime();
    }

    // Calculate minimum and maximum
    if (timer.cpuTime() < minimum)
      minimum = timer.cpuTime();
    if (timer.cpuTime() > maximum)
      maximum = timer.cpuTime();

    // Write the output
    cout << counter << " inversions using inverse() took place\n";
    cout << "Time is " << timer.cpuTime() << " secs\n";
    cout << "Inversion time is " << 1000000 * timer.cpuTime() / N << " µsecs\n";

    // Reset the timer and the counter
    timer.reset();
    counter = 0;
  }
}

// Write the minimum and maximum values
void Summary(const double& minimum, const double& maximum)
{
  cout << "\nThe minimum time for 1 run was " << minimum << " secs\n";
  cout << "The maximum time for 1 run was " << maximum << " secs\n";
}


int main()
{
  // What package are we using?
  Identification();

  // Get the dim(-ension)
  int dim = 0; GetDimension(dim);

  // Define a "dim x dim" matrix
  Matrix m (dim, dim);  

  // Define how the matrix was declared
  Describe (dim);

  // Build the "dim x dim" matrix
  BuildMatrix(m, dim);

  // Declare the matrix specialized?
  DeclareMatrixSpecialized(m, dim);

  // Show the actual matrix
  Show(m);

  // Inversion parameters
  int inv = 0; GetCalcs(inv); // number of inversion calculations
  int run = 0; GetRuns(run);  // number of runs
  double min = 0.;  // minimum value
  double max = 0.;  // maxinum value

  // Invert the matrix
  Calculate(m, inv, run, min, max);

  // Display the minimum and maximum
  Summary(min, max);

  return 0;
}