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

// CLHEP: HepSymMatrix inversion

#include "CLHEP/Matrix/SymMatrix.h"

using namespace std;

// Function declarations
void Identification();
void GetDimension(int&);
void Describe(const int&);
void BuildMatrix(HepSymMatrix&, const int&);
void Show(const HepSymMatrix&);
void GetMethod(char&);
void GetCalcs(int&);
void GetRuns(int&);
void Calculate(HepSymMatrix&, const int&, const char&, const int&, const int&);


// Package ID
void Identification()
{
  cout << "Using CLHEP\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 'HepSymMatrix (" << dimension << ")'\n";
}

// Build the matrix
void BuildMatrix(HepSymMatrix& matrix, const int& dimension)
{
  for (int i=0; i<dimension; ++i)
  {
    for (int j=0; j<dimension; ++j)
    {
      if (i == j)
        matrix[i][j] = 1000 + i*100;
      else if (i < j)
        matrix[i][j] = 2 + i*(dimension*2) + j*2;
      else if (j < i)
        matrix[i][j] = 2 + i*2 + j*(dimension*2);
    }
  }
}

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

// Get the inversion method
void GetMethod(char& method)
{
  cout << "Which inversion method you want to use?\n";
  cout << "Enter s for inverSe(), t for inverT(): ";
  cin >> method;

  while ( (method != 's') && (method != 'S') &&
          (method != 't') && (method != 'T') )
  {
    cout << "Enter s or t: ";
    cin >> method;
  }
}

// Get the number of calculations
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(HepSymMatrix& matrix, const int& dimension, const char& imethod,
               const int& N, const int& runs)
{
  // Define a counter
  int counter = 0;

  // Define a matrix to hold the result
  HepSymMatrix inverse (dimension);

  // Do we have an SPD matrix?
  int ifail = 0;

  if ( (imethod == 's') || (imethod == 'S') )
  {
    // Loop "run" times
    for (int r=0; r<runs; ++r)
    {
      // Do the inversions
      for (int i=0; i<N; ++i)
      {
        // The contents of the matrix are not replaced by the inverse
	inverse = matrix.inverse(ifail);
	// thus this is necessary
	matrix = inverse;
        ++counter;
      }

      // Write the output
      cout << counter << " inversions using inverse() took place\n";
      cout << "Inverse =" << inverse << "\n";

      // Reset the counter
      counter = 0;
    }
  }
  else if ( (imethod == 't') || (imethod == 'T') )
  {
    // Loop "run" times
    for (int r=0; r<runs; ++r)
    {
      // Do the inversions
      for (int i=0; i<N; ++i)       
      {
        // The contents of the matrix are replaced by the inverse
	matrix.invert(ifail);
	// thus this isn't necessary
	// matrix = inverse;
        ++counter;
      }

      // Write the output
      cout << counter << " inversions using inverse() took place\n";
      cout << "Inverse =" << matrix << "\n";

      // Reset the counter
      counter = 0;
    }
  }
}


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

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

  // Define a "dim x dim" HepSymMatrix
  HepSymMatrix m (dim);

  // Describe how the matrix was defined
  Describe (dim);

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

  // Show the actual matrix
  Show(m);

  // Inversion parameters
  char meth = ' '; GetMethod(meth); // inversion method
  int tim = 0; GetCalcs(tim);       // number of inversions
  int run = 0; GetRuns(run);        // number of runs

  // Invert the matrix
  Calculate(m, dim, meth, tim, run);

  return 0;
}