#include <iostream>
#include <complex>
#include "vector.h"
#include "matrix.h"
#include "timer.h"

// another frend operator to compute eigenvalues
inline int Eigenvalues(Matrix& H, Vector& Er, Vector& Ei)
{
  static Vector work(4*H.size0());
  if (H.size0()!=H.size1()) {std::cerr<<"Can not diagonalize non-square matrix!"<<std::endl;}
  work.resize(4*H.size0());
  Er.resize(H.size0());
  Ei.resize(H.size0());
  xgeev_(H.size0(), H.m, H.size1(), &(Er[0]), &(Ei[0]), &(work[0]), work.size());
}

Vector operator*(const Matrix& A, const Vector& b)
{
  assert(A.size1()==b.size()); // sizes match
  Vector c(A.size0());
  for (int i=0; i<A.size0(); i++){
    double sum=0;
    for (int j=0; j<A.size1(); j++) sum += A(i,j)*b[j];
    c[i] = sum;
  }
  return c;
}

using namespace std;

void TestMatrixVector(int N, int M)
{
  cout<<"*************************************************"<<endl;
  cout<<"Test Matrix - Vector operation"<<endl;
  
  Matrix A(N,M);
  A.random();

  cout<<"Matrix A is "<<endl<<A;
  
  Vector b(M);
  b.random();

  cout<<"Vector b is "<<b<<endl;
  
  Vector c = A*b;

  cout<<"Product A*b is "<<c<<endl;
  cout<<endl;
}

void TestMatrixMatrix(int N0)
{
  cout<<"*************************************************"<<endl;
  cout<<"Test Matrix - Matrix operation"<<endl;
  
  Matrix A(N0,N0), B(N0,N0);

  A.random();
  B.random();
  
  cout<<"Matrix A is "<<endl;
  cout<<A<<endl;
  cout<<"Matrix B is "<<endl;
  cout<<B<<endl;
  
  Matrix C;
  //  C.Product(A,B);
  C = A * B;
  cout<<"Matrix A*B is "<<endl<<C;


  Vector Er(A.size0()), Ei(A.size0());
  
  Eigenvalues(A, Er, Ei);
  cout<<"Eigevalues of A are: ";
  for (int i=0; i<Er.size(); i++) cout<<Er[i]<<" "<<Ei[i]<<"*i, ";
  cout<<endl;

  complex<double> detA = 1;
  for (int i=0; i<Er.size(); i++) detA *= complex<double>(Er[i],Ei[i]);
    
  Eigenvalues(B, Er, Ei);
  cout<<"Eigevalues of B are: ";
  for (int i=0; i<Er.size(); i++) cout<<Er[i]<<" "<<Ei[i]<<"*i, ";
  cout<<endl;

  complex<double> detB = 1;
  for (int i=0; i<Er.size(); i++) detB *= complex<double>(Er[i],Ei[i]);
  
  Eigenvalues(C, Er, Ei);
  cout<<"Eigevalues of C are: ";
  for (int i=0; i<Er.size(); i++) cout<<Er[i]<<" "<<Ei[i]<<"*i, ";
  cout<<endl;
  
  complex<double> detC = 1;
  for (int i=0; i<Er.size(); i++) detC *= complex<double>(Er[i],Ei[i]);

  cout<<"Determinant of product is "<<detC<<" while product of determinants is "<<detA*detB<<endl;
  cout<<endl;
}

void TestSpeed()
{
  cout<<"*************************************************"<<endl;
  cout<<" Test speed"<<endl;
  int N0=2;
  while(N0<2000){
    N0 = static_cast<int>(N0*1.5);
      
    Matrix A(N0,N0), B(N0,N0);
    Matrix C;

    A.random();
    B.random();
    
    int Nmax = static_cast<int>(1e8/(N0*N0*N0));
    if (Nmax>1e6) Nmax = static_cast<int>(1e6);
    if (Nmax<1) Nmax = 1;
    
    Timer t;
    for (int i=0; i<Nmax; i++) C = A * B; //C.Product(A,B);//C=A*B;//
    t.stop();
  
    double time = t.toSeconds()/Nmax;
    double time_per_FPO = time/(2.*N0*N0*N0);
    double Mflops = 1e-6/time_per_FPO;
  
    cout<<"It takes "<<setw(14)<<time<<"s to multiply matrix of size "<<setw(4)<<N0<<" efficiency is "<<setw(5)<<static_cast<int>(Mflops)<<" Mflops"<<endl;
    //  compare with cat /proc/cpuinfo
  }
}

int main()
{
  srand48(time(0));
  
  TestMatrixVector(5, 6);
  
  TestMatrixMatrix(5);
    
  TestSpeed();
  
  return 0;
}