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

// CLHEP: HepVector*HepMatrix multiplication
// (In order to perform vector*matrix in CLHEP we have to transpose
// the vector. However, when we transpose a HepVector we get back a
// HepMatrix. Thus, we end up performing matrix*matrix multiplication
// and not vector*matrix.)

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

using namespace std;

// Function declarations
void Identification();
void GetDimensions(int&);
void Describe(const int&);
void BuildMatrix(HepMatrix&, const int&);
void BuildVector(HepVector&, const int&);
void Show(const HepMatrix&, const HepVector&);
void GetCalcs(int&);
void GetRuns(int&);
void Calculate(HepMatrix&, HepVector&, const int&, const int&, const int&);


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

// Ask for the desired dimensions
void GetDimensions(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 and the vector were defined
void Describe(const int& dimension)
{
  cout << "Matrix was defined 'HepMatrix ("
       << dimension << ", " << dimension << ")'\n";
  cout << "Vector was defined 'HepVector (" << dimension << ")'\n";
}

// Build the matrix
void BuildMatrix(HepMatrix& 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 (i > j)
        matrix[i][j] = 2 + i*2 + j*(dimension*2);
    }
  }
}

// Build the vector
void BuildVector(HepVector& vector, const int& dimension)
{
  for (int i=0; i<dimension; ++i)
  {
    vector[i] = 2*(i+1);
  }
}

// Show the matrix and the vector
void Show(const HepMatrix& matrix, const HepVector& vector)
{
  cout << "Using matrix =" << matrix << "\n";
  cout << "Using vector =" << vector.T() << "\n";     
  // doing .T() here doesn't change the vector
}

// 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 HepMatrix with the HepVector
void Calculate(HepMatrix& matrix, HepVector& vector, const int& dimension,
               const int& N, const int& runs)
{  
  // Define a counter
  int counter = 0;

  // Define a HepMatrix to hold the transpose of the HepVector
  HepMatrix transpose (dimension, dimension);

  // Transpose the HepVector so you can do vector*matrix
  transpose = vector.T();

  // Define a vector to hold the result
  HepVector product_vector (dimension);

  // Define a HepMatrix to hold the transpose of the product_vector
  HepMatrix product (dimension, dimension);

  // Transpose the product vector
  product = product_vector.T();

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

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

    // Reset the counter
    counter = 0;
  }
}


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

  // Get the dim(-ensions)
  int dim = 0; GetDimensions(dim);

  // Define a "dim x dim" HepMatrix
  HepMatrix m (dim, dim);
  
  // Define a "dim-element" HepVector
  HepVector v (dim);

  // Describe how the matrix and the vector were defined
  Describe (dim);

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

  // Build the "dim-element" vector
  BuildVector(v, dim);

  // Show the actual matrix/vector
  Show(m, v);

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

  // Multiply the matrix with the vector
  Calculate(m, v, dim, calc, run);

  return 0;
}