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

// CLHEP: HepVector*HepVector multiplication
//        vector form like in RowVector*ColumnVector
// (In order to perform vector*vector in CLHEP we have to transpose
// one of the two vectors. Here we transpose the first vector.
// However, when we transpose a HepVector we get back a HepMatrix.
// Thus, we end up performing matrix*vector and not vector*vector.)

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

using namespace std;

// Function declarations
void Identification();
void GetDimension(int&);
void Describe(const int&);
void BuildVectors(HepVector&, HepVector&, const int&);
void Show(const HepVector&, const HepVector&);
void GetCalcs(int&);
void GetRuns(int&);
void Calculate(HepVector&, HepVector&, 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 two vectors were defined
void Describe(const int& dimension)
{
  cout << "Vectors were defined 'HepVector (" << dimension << ")'\n";
}

// Build the two vectors
void BuildVectors(HepVector& vector1, HepVector& vector2, const int& dimension)
{
  for (int i=0; i<dimension; ++i)
  {
    vector1[i] = 2*(i+1);
    vector2[i] = 2*(i+1);
  }
} 

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

// Get the number of multiplications
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 two vectors
void Calculate(HepVector& vector1, HepVector& vector2,
               const int& N, const int& runs)
{
  // Define a counter
  int counter = 0;

  // Define a HepMatrix to hold the inverted vector
  HepMatrix vector1_inverted;

  // Invert the second vector
  vector1_inverted = vector1.T();

  // Define a vector to hold the result
  HepVector product(1);

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Do the multiplications
    for (int i=0; i<N; ++i)
    {
      product += vector1_inverted * vector2;
      ++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(-ension)
  int dim = 0; GetDimension(dim);
  
  // Definee two "dim-element" HepVectors
  HepVector v1 (dim);
  HepVector v2 (dim);

  // Describe how the two vectors were declared
  Describe (dim);

  // Build the vectors
  BuildVectors(v1, v2, dim);

  // Show the actual vectors
  Show(v1, v2);

  // Multiplication parameters
  int mul = 0; GetCalcs(mul); // number of multiplications
  int run = 0; GetRuns(run);  // number of runs

  // Multiply the two vectors
  Calculate(v1, v2, mul, run);

  return 0;
}