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

// CovMatrices: CovMatrixX*CovMatrixX multiplication

#include "CovMatrices/CovMatrix.h"
#include "CovMatrices/CovMatrix2.h"
#include "CovMatrices/CovMatrix3.h"
#include "CovMatrices/CovMatrix4.h"
#include "CovMatrices/CovMatrix5.h"
#include "CovMatrices/CovMatrix6.h"
#include <vector>
#include "CovMatrices/Vector2.h"
#include "CovMatrices/Vector3.h"
#include "CovMatrices/Vector4.h"
#include "CovMatrices/Vector5.h"
#include "CovMatrices/Vector6.h"

using namespace std;

using CovMatrices::CovMatrix;
using CovMatrices::CovMatrix2;
using CovMatrices::CovMatrix3;
using CovMatrices::CovMatrix4;
using CovMatrices::CovMatrix5;
using CovMatrices::CovMatrix6;
using CovMatrices::Vector2;
using CovMatrices::Vector3;
using CovMatrices::Vector4;
using CovMatrices::Vector5;
using CovMatrices::Vector6;

// Function declarations
void Identification();
void GetDimension(int&);
void Describe(const int&);
void BuildCovMatrix(CovMatrix&);
void BuildCovMatrix2(CovMatrix2&);
void BuildCovMatrix3(CovMatrix3&);
void BuildCovMatrix4(CovMatrix4&);
void BuildCovMatrix5(CovMatrix5&);
void BuildCovMatrix6(CovMatrix6&);
void Show(const CovMatrix&, const CovMatrix&);
void Show2(const CovMatrix2&, const CovMatrix2&);
void Show3(const CovMatrix3&, const CovMatrix3&);
void Show4(const CovMatrix4&, const CovMatrix4&);
void Show5(const CovMatrix5&, const CovMatrix5&);
void Show6(const CovMatrix6&, const CovMatrix6&);
int GetCalcs();
int GetRuns();
void Calculate(CovMatrix&, CovMatrix&, const int&, const int&);
void Calculate2(CovMatrix2&, CovMatrix2&, const int&, const int&);
void Calculate3(CovMatrix3&, CovMatrix3&, const int&, const int&);
void Calculate4(CovMatrix4&, CovMatrix4&, const int&, const int&);
void Calculate5(CovMatrix5&, CovMatrix5&, const int&, const int&);
void Calculate6(CovMatrix6&, CovMatrix6&, const int&, const int&);


// Package ID
void Identification()
{
  cout << "Using CovMatrices\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)
{
  if (dimension == 1)
    cout << "Matrices were defined 'CovMatrix (" << dimension << ")'\n";
  else
    cout << "Matrices defined 'CovMatrix" << dimension << "'\n";
}

// Built a 1-dimension CovMatrix
void BuildCovMatrix(CovMatrix& matrix)
{
  matrix(0,0) = 1000;
}

// Built a 2-dimension CovMatrix
void BuildCovMatrix2(CovMatrix2& matrix2)
{
  for (int i=0; i<2; ++i)
  {
    for (int j=0; j<2; ++j)
    {
      if (i == j)
        matrix2(i,j) = 1000 + i*100;
      else if (i > j)
        matrix2(i,j) = 2 + i*2 + j*(2*2);
    }
  }
}

// Built a 3-dimension CovMatrix
void BuildCovMatrix3(CovMatrix3& matrix3)
{
  for (int i=0; i<3; ++i)
  {
    for (int j=0; j<3; ++j)
    {
      if (i == j)
        matrix3(i,j) = 1000 + i*100;
      else if (i > j)
        matrix3(i,j) = 2 + i*2 + j*(3*2);
    }
  }
}

// Built a 4-dimension CovMatrix
void BuildCovMatrix4(CovMatrix4& matrix4)
{
  for (int i=0; i<4; ++i)
  {
    for (int j=0; j<4; ++j)
    {
      if (i == j)
        matrix4(i,j) = 1000 + i*100;
      else if (i > j)
        matrix4(i,j) = 2 + i*2 + j*(4*2);
    }
  }
}

// Built a 5-dimension CovMatrix
void BuildCovMatrix5(CovMatrix5& matrix5)
{
  for (int i=0; i<5; ++i)
  {
    for (int j=0; j<5; ++j)
    {
      if (i == j)
        matrix5(i,j) = 1000 + i*100;
      else if (i > j)
        matrix5(i,j) = 2 + i*2 + j*(5*2);
    }
  }
}

// Built a 6-dimension CovMatrix
void BuildCovMatrix6(CovMatrix6& matrix6)
{
  for (int i=0; i<6; ++i)
  {
    for (int j=0; j<6; ++j)
    {
      if (i == j)
        matrix6(i,j) = 1000 + i*100;
      else if (i > j)
        matrix6(i,j) = 2 + i*2 + j*(6*2);
    }
  }
}

// Show the 1-dimensional CovMatrix
void Show(const CovMatrix& matrix1, const CovMatrix& matrix2)
{
  cout << "Using matrix =\n" << matrix1 << "\n";
  cout << "Using matrix =\n" << matrix2 << "\n";
}

// Show the 2-dimensional CovMatrix
void Show2(const CovMatrix2& matrix1, const CovMatrix2& matrix2)
{
  cout << "Using matrix =\n" << matrix1 << "\n";
  cout << "Using matrix =\n" << matrix2 << "\n";
}

// Show the 3-dimensional CovMatrix
void Show3(const CovMatrix3& matrix1, const CovMatrix3& matrix2)
{
  cout << "Using matrix =\n" << matrix1 << "\n";
  cout << "Using matrix =\n" << matrix2 << "\n";
}

// Show the 4-dimensional CovMatrix
void Show4(const CovMatrix4& matrix1, const CovMatrix4& matrix2)
{
  cout << "Using matrix =\n" << matrix1 << "\n";
  cout << "Using matrix =\n" << matrix2 << "\n";
}

// Show the 5-dimensional CovMatrix
void Show5(const CovMatrix5& matrix1, const CovMatrix5& matrix2)
{
  cout << "Using matrix =\n" << matrix1 << "\n";
  cout << "Using matrix =\n" << matrix2 << "\n";
}

// Show the 6-dimensional CovMatrix
void Show6(const CovMatrix6& matrix1, const CovMatrix6& matrix2)
{
  cout << "Using matrix =\n" << matrix1 << "\n";
  cout << "Using matrix =\n" << matrix2 << "\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 1-dimensional matrices
void Calculate(CovMatrix& matrix1, CovMatrix& matrix2,
               const int& N, const int& runs)
{
  // Define a counter
  int counter = 0;

  // Define a matrix to hold the result
  CovMatrix product(1);

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Do the multiplications ("matrix version" of vector dot product)
    for (int n=0; n<N; ++n)
    {
      for (int i=0; i<1; ++i)
      {
        for (int j=0; j<1; ++j)
	{
	  product(i,j) =    matrix1(i,0) * matrix2(j,0);
	}
      }
      ++counter;
      matrix1 = product;
    }

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

    // Reset the counter
    counter = 0;
  }
}

// Multiply the 2-dimensional matrix and vector
void Calculate2(CovMatrix2& matrix1, CovMatrix2& matrix2,
                const int& N, const int& runs)
{
  // Define a counter
  int counter2 = 0;

  // Define a matrix to hold the result
  CovMatrix2 product2;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Do the multiplications ("matrix version" of vector dot product)
    for (int n=0; n<N; ++n)
    {
      for (int i=0; i<2; ++i)
      {
        for (int j=0; j<2; ++j)
	{
	  product2(i,j) =   matrix1(i,0) * matrix2(j,0)
	                  + matrix1(i,1) * matrix2(j,1);
	}
      }
      ++counter2;
      matrix1 = product2;
    }

    // Write the output
    cout << counter2 << " multiplications took place\n";
    cout << "Product =\n" << product2 << "\n";

    // Reset the counter
    counter2 = 0;
  }
}
 
// Multiply the 3-dimensional matrices
void Calculate3(CovMatrix3& matrix1, CovMatrix3& matrix2,
                const int& N, const int& runs)
{
  // Define a counter
  int counter3 = 0;

  // Define a matrix to hold the result
  CovMatrix3 product3;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Do the multiplications ("matrix version" of vector dot product)
    for (int n=0; n<N; ++n)
    {
      for (int i=0; i<3; ++i)
      {
        for (int j=0; j<3; ++j)
	{
	  product3(i,j) =   matrix1(i,0) * matrix2(j,0)
	                  + matrix1(i,1) * matrix2(j,1)
	                  + matrix1(i,2) * matrix2(j,2);
	}
      }
      ++counter3;
      matrix1 = product3;
    }

    // Write the output
    cout << counter3 << " multiplications took place\n";
    cout << "Product =\n" << product3 << "\n";

    // Reset the counter
    counter3 = 0;
  }
}

// Multiply the 4-dimensional matrices
void Calculate4(CovMatrix4& matrix1, CovMatrix4& matrix2,
                const int& N, const int& runs)
{
  // Define a counter
  int counter4 = 0;

  // Define a matrix to hold the result
  CovMatrix4 product4;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Do the multiplications ("matrix version" of vector dot product)
    for (int n=0; n<N; ++n)
    {
      for (int i=0; i<4; ++i)
      {
        for (int j=0; j<4; ++j)
	{
	  product4(i,j) =   matrix1(i,0) * matrix2(j,0)
	                  + matrix1(i,1) * matrix2(j,1)
	                  + matrix1(i,2) * matrix2(j,2)
	                  + matrix1(i,3) * matrix2(j,3);
	}
      }
      ++counter4;
      matrix1 = product4;
    }

    // Write the output
    cout << counter4 << " multiplications took place\n";
    cout << "Product =\n" << product4 << "\n";

    // Reset the counter
    counter4 = 0;  
  }
}

// Multiply the 5-dimensional matrices
void Calculate5(CovMatrix5& matrix1, CovMatrix5& matrix2,
                const int& N, const int& runs)
{
  // Define a counter
  int counter5 = 0;

  // Define a matrix to hold the result
  CovMatrix5 product5;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Do the multiplications ("matrix version" of vector dot product)
    for (int n=0; n<N; ++n)
    {
      for (int i=0; i<5; ++i)
      {
        for (int j=0; j<5; ++j)
	{
	  product5(i,j) =   matrix1(i,0) * matrix2(j,0)
	                  + matrix1(i,1) * matrix2(j,1)
	                  + matrix1(i,2) * matrix2(j,2)
	                  + matrix1(i,3) * matrix2(j,3)
	                  + matrix1(i,4) * matrix2(j,4);
	}
      }
      ++counter5;
      matrix1 = product5;
    }

    // Write the output
    cout << counter5 << " multiplications took place\n";
    cout << "Product =\n" << product5 << "\n";

    // Reset the counter
    counter5 = 0;
  }
}

// Multiply the 6-dimensional matrices
void Calculate6(CovMatrix6& matrix1, CovMatrix6& matrix2,
                const int& N, const int& runs)
{
  // Define a counter
  int counter6 = 0;

  // Define a matrix to hold the result
  CovMatrix6 product6;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Do the multiplications ("matrix version" of vector dot product)
    for (int n=0; n<N; ++n)
    {
      for (int i=0; i<6; ++i)
      {
        for (int j=0; j<6; ++j)
	{
	  product6(i,j) =   matrix1(i,0) * matrix2(j,0)
	                  + matrix1(i,1) * matrix2(j,1)
	                  + matrix1(i,2) * matrix2(j,2)
	                  + matrix1(i,3) * matrix2(j,3)
	                  + matrix1(i,4) * matrix2(j,4)
			  + matrix1(i,5) * matrix2(j,5);
          ++counter6;
	}
      }
      ++counter6;
      matrix1 = product6;
    }

    // Write the output
    cout << counter6 << " multiplications took place\n";
    cout << "Product =\n" << product6 << "\n";

    // Reset the counter
    counter6 = 0;
  }
}


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

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

  // Define the CovMatrices
  CovMatrix m11(1), m12(1);
  CovMatrix2 m21, m22;
  CovMatrix3 m31, m32;
  CovMatrix4 m41, m42;
  CovMatrix5 m51, m52;
  CovMatrix6 m61, m62;

  // Describe how the matrices were defined
  Describe (dim);

  // Build the matrices
  switch(dim)
  {
    case 1:
      BuildCovMatrix(m11);
      BuildCovMatrix(m12);
      break;
    case 2:
      BuildCovMatrix2(m21);
      BuildCovMatrix2(m22);
      break;
    case 3:
      BuildCovMatrix3(m31);
      BuildCovMatrix3(m32);
      break;
    case 4:
      BuildCovMatrix4(m41);
      BuildCovMatrix4(m42);
      break;
    case 5:
      BuildCovMatrix5(m51);
      BuildCovMatrix5(m52);
      break;
    case 6:
      BuildCovMatrix6(m61);
      BuildCovMatrix6(m62);
      break;
  }

  // Show the actual matrices
  switch(dim)
  {
    case 1:
      Show(m11, m12);
      break;
    case 2:
      Show2(m21, m22);
      break;
    case 3:
      Show3(m31, m32);
      break;
    case 4:
      Show4(m41, m42);
      break;
    case 5:
      Show5(m51, m52);
      break;
    case 6:
      Show6(m61, m62);
      break;
  }

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

  // Do the multiplications
  switch(dim)
  {
    case 1:
      Calculate(m11, m12, calc, run);
      break;
    case 2:
      Calculate2(m21, m22, calc, run);
      break;
    case 3:
      Calculate3(m31, m32, calc, run);
      break;
    case 4:
      Calculate4(m41, m42, calc, run);
      break;
    case 5:
      Calculate5(m51, m52, calc, run);
      break;
    case 6:
      Calculate6(m61, m62, calc, run);
      break;
  }

  return 0;
}