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

// CLHEP: HepSymMatrix multiplication

#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"

using namespace std;

// Function declarations
void GetDimension();
void Describe(const int&);
void BuildMatrix(HepSymMatrix&, const int&);
void Show(const HepSymMatrix&);
void GetCalcs(int&);
void GetRuns(int&);
void Calculate(HepSymMatrix&, HepSymMatrix&,
               const int&, 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 matrices were defined
void Describe(const int& dimension)
{
  cout << "Matrices were 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 number of calculations
void GetCalcs(int& N)
{
  cout << "How many multiplications? ";
  cin >> N;
}

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

// Multiply the matrices
void Calculate(HepSymMatrix& matrix1, HepSymMatrix& matrix2,
               const int& dimension, const int& N, const int& runs)
{
  // Define a counter
  int counter = 0;

  // Define a matrix to hold the result
  // (define it HepMatrix and not HepSymMatrix)
  HepMatrix product (dimension, dimension); 

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  { 
    // Do the multiplications
    for (int i=0; i<N; ++i)
    {
      product = matrix1 * matrix2;

      // Since a HepMatrix (product) can't be assigned to
      // a HepSymMatrix (matrix1) we assign every
      // individual element from the one matrix to the other
      for(int x=0; x<dimension; ++x)
      {
        for(int y=0; y<dimension; ++y)
	{
	  matrix1[x][y] = product[x][y];
	}
      }      
      ++counter;
    }

    // Write the output
    cout << counter << " multiplications took place\n";
    cout << "Product =" << product;

    // Reset the counter
    counter = 0;
  }
}


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

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

  // Define two "dim x dim" matrices
  HepSymMatrix m1 (dim);
  HepSymMatrix m2 (dim);

  // Define how the matrices were declared
  Describe (dim);

  // Build the matrices
  BuildMatrix(m1, dim);
  BuildMatrix(m2, dim);

  // Show the actual matrices
  Show(m1);
  Show(m2);

  // Calculation parameters
  int tim = 0; GetCalcs(tim); // times to perform the calculation
  int run = 0; GetRuns(run);  // number of runs

  // Do the matrix multiplications
  Calculate(m1, m2, dim, tim, run);

  return 0;
}