// comparison of matrix products 
//  armadillo vs. plain loop
// Optimize as well as you can, e.g.
//g++ -O3 -funroll-all-loops arma_matrix_speed2.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=500,m1=500,m2=100;
  // armadillo
  mat A = randu<mat>(n,m1); // mat is double
  mat B = randu<mat>(m1,m2);  
  mat C = randu<mat>(m2,n);  
  mat D(n,n);
  tic=clock();
  D = A*B*C;
  toc=clock();
  cout <<"       armadillo timing = "<< (toc-tic)/time_unit<<" sec"<<endl;

  /*
  // plain for loop using armadillo matrices
  tic=clock();
  for (unsigned i = 0; i<n; ++i)
    for (unsigned j = 0; j<n; ++j)
      for (unsigned k = 0; k<m1; ++k)
	for (unsigned l = 0; l<m2; ++l)
	  D(i,j) = A(i,k)*B(k,l)*C(l,j);	
  toc=clock();
  cout <<"  armadillo loop timing = "<< (toc-tic)/time_unit<<" sec"<<endl;
  */

  // ordinary data  version
  double AA[n][m1],BB[m1][m2],CC[m2][n],DD[n][n];
  double rand_unit=static_cast<double>(RAND_MAX);
  for (unsigned i = 0; i<n; ++i) 
    for (unsigned j = 0; j<m1; ++j) 
      AA[i][j] = A(i,j) ;// rand()/rand_unit;
  for (unsigned i = 0; i<m1; ++i) 
    for (unsigned j = 0; j<m2; ++j) 
      BB[i][j] = B(i,j) ;// rand()/rand_unit;
  for (unsigned i = 0; i<m2; ++i) 
    for (unsigned j = 0; j<n; ++j) 
      CC[i][j] = C(i,j);//rand()/rand_unit;
  for (unsigned i = 0; i<n; ++i) 
    for (unsigned j = 0; j<n; ++j) 
      DD[i][j]=0;

  tic=clock();
  for (unsigned i = 0; i<n; ++i){
    for (unsigned j = 0; j<n; ++j){
      for (unsigned k = 0; k<m1; ++k){
	for (unsigned l= 0; l<m2; ++l){
	  DD[i][j] += AA[i][k]*BB[k][l]*CC[l][j];	
	}
      }
    }
  }
  toc=clock();
  cout<<endl;
  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<<D(i,i)<<" ";}
  cout<<endl;
  for (unsigned i=0; i<10;++i) {cout<<DD[i][i]<<" ";}
  cout<<endl;
  return 0;
}