// comparison of matrix products 
//  armadillo vs. plain loop
// Optimize as well as you can, e.g.
// g++ -O3 -funroll-all-loops arma_matrix_speed.cpp -L/usr/lib64/atlas -lcblas
#include <iostream>
#include <armadillo>
#include <ctime>

using namespace std;
using namespace arma;

int main()
{
  clock_t tic,toc;
  double time_unit=static_cast<double>(CLOCKS_PER_SEC);
  unsigned n=1000,m=1200;
  // armadillo
  mat A = randu<mat>(n,m); // mat is double
  mat B = randu<mat>(m,n);  
  mat C(n,n);
  tic=clock();
  C = A*B;
  toc=clock();
  cout <<"       armadillo timing = "<< (toc-tic)/time_unit<<" sec"<<endl;
 
  // plain for loop using armadillo matrices
  // extremely slow!
  tic=clock();
  C.zeros(n,n);
  for (unsigned i = 0; i<n; ++i)
    for (unsigned j = 0; j<n; ++j)
      for (unsigned k = 0; k<m; ++k)
	C(i,j) += A(i,k)*B(k,j);	  
  // try: C(i,j) = C(i,j) + A(i,k)*B(k,j);
  toc=clock();
  cout <<"  armadillo loop timing = "<< (toc-tic)/time_unit<<" sec"<<endl;
 

  // ordinary data  version
  double AA[n][m],BB[m][n],CC[n][n];
  for (unsigned i = 0; i<n; ++i) {
    for (unsigned j = 0; j<m; ++j) {
      AA[i][j] = A(i,j);
      BB[j][i] = B(j,i);
    }
  }
  for (unsigned i = 0; i<n; ++i) 
    for (unsigned j = 0; j<n; ++j) 
      CC[i][j]=0;

  tic=clock();
  for (unsigned i = 0; i<n; ++i){
    for (unsigned j = 0; j<n; ++j){
      for (unsigned k = 0; k<m; ++k)
	CC[i][j] += AA[i][k]*BB[k][j];	 // this += makes this one slow 
      //CC[i][j] = AA[i][k]*BB[k][j];  // this is rather fast
    }
  }
  toc=clock();
  cout <<"      plain loop timing = "<< (toc-tic)/time_unit<<" sec"<<endl;
  cout<<endl;
  cout<<"a few matrix values to check results"<<endl;
  for (unsigned i=0; i<10;++i) {cout<<C(i,i)<<" ";}
  cout<<endl;
  for (unsigned i=0; i<10;++i) {cout<<CC[i][i]<<" ";}
  cout<<endl;



  return 0;
}