#include <time.h>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/fwd.hpp>
#include <boost/numeric/ublas/expression_types.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/io.hpp>

int main()
{
  using namespace boost::numeric::ublas;
  int i,j,ta1,ta2,dt_normal,tb1,tb2,dt_Sparse,tc1,tc2,dt_Symmetry,dt_smc,dt_smr;
  int dt_1,dt_2,dt_3,dt_4,dt_5,dt_6,dt_7,dt_8,dt_9,dt_10,dt_11;

  const int size = 800;
  const int count = 2;

  {
    symmetric_matrix<float,upper> sy1(size,size);
    symmetric_matrix<float,lower> sy2(size, size),sy3(size,size);
    compressed_matrix<float,column_major> sp1(size,size),sp2(size,size),sp3(size,size); 
    matrix <float, column_major> m4(size,size);
    matrix <float, row_major> m5(size,size),m6(size,size);

    //---initial matrix-------//
    for(i=0;i<size;i++)
      {
        for(j=0;j<size-1;j++)
          {               
            //           sp1(i,j)=0;
            //           sp2(i,j)=0;
            sy1(i,j)=0;
            sy2(i,j)=0;
            m4(i,j)=0;
            m5(i,j)=0;
          }
        sp1(i,j)=1;
        sp2(i,j)=1;
        sy1(i,j)=1;
        sy2(i,j)=1;
        m4(i,j)=1;
        m5(i,j)=1;
      }

    //--sparse * sparse column mode //...(1)
    tb1=clock();
    for(i=0;i<count;i++)
      {
        sp3=prod(sp1,sp2);
      }
    tb2=clock();
    dt_1=tb2-tb1;   
    printf("%d\n", dt_1);

    //--sparse * normal column mode //...(2)
    ta1=clock();
    for(i=0;i<count;i++)
      {
        sp3=prod(sp1,m4);
      }
    ta2=clock();
    dt_2=ta2-ta1;
    printf("%d\n", dt_2);
    //--sparse * normal row mode //...(3)
    ta1=clock();
    for(i=0;i<count;i++)
      {
        sp3=prod(sp1,m5);
      }
    ta2=clock();
    dt_3=ta2-ta1;
    printf("%d\n", dt_3);

    //--sparse * symmerty upper mode //...(4)
    ta1=clock();
    for(i=0;i<count;i++)
      {
        sp3=prod(sp1,sy1);
      }
    ta2=clock();
    dt_4=ta2-ta1;
    printf("%d\n", dt_4);

    //--sparse * symmerty lower mode //...(5)
    ta1=clock();
    for(i=0;i<count;i++)
      {
        sp3=prod(sp1,sy2);
      }
    ta2=clock();
    dt_5=ta2-ta1;
    printf("%d\n", dt_5);

    //--dense * dense  //...(5)
    {
      int ta1=clock();
      for(i=0;i<count;i++)
        {
          m6=prod(m4,m5);
        }
      int ta2=clock();
      dt_5=ta2-ta1;
      printf("%d\n", dt_5);
    }
    //--sparse * sparse column mode //...(1)
    if (0) {
      tb1=clock();
      for(i=0;i<count;i++)
        {
          m6=prod(sp1,sp2);
        }
      tb2=clock();
      dt_1=tb2-tb1;   
      printf("%d\n", dt_1);
    }
  }

  {
    symmetric_matrix<float,upper> sy1(size,size);
    symmetric_matrix<float,lower> sy2(size, size),sy3(size,size);
    compressed_matrix<float,column_major> sp1(size,size),sp2(size,size),sp3(size,size); 
    matrix <float, column_major> m4(size,size);
    matrix <float, row_major> m5(size,size),m6(size,size);

    //---initial matrix-------//
    for(i=0;i<size;i++)
      {
        for(j=0;j<size-1;j++)
          {               
            //           sp1(i,j)=0;
            //           sp2(i,j)=0;
            sy1(i,j)=0;
            sy2(i,j)=0;
            m4(i,j)=0;
            m5(i,j)=0;
          }
        sp1(i,j)=1;
        sp2(i,j)=1;
        sy1(i,j)=1;
        sy2(i,j)=1;
        m4(i,j)=1;
        m5(i,j)=1;
      }

    //--sparse * sparse column mode //...(1)
    tb1=clock();
    for(i=0;i<count;i++)
      {
        axpy_prod(sp1,sp2,sp3,true);
      }
    tb2=clock();
    dt_1=tb2-tb1;   
    printf("%d\n", dt_1);

    //--sparse * normal column mode //...(2)
    ta1=clock();
    for(i=0;i<count;i++)
      {
        axpy_prod(sp1,m4,sp3,true);
      }
    ta2=clock();
    dt_2=ta2-ta1;
    printf("%d\n", dt_2);
    //--sparse * normal row mode //...(3)
    ta1=clock();
    for(i=0;i<count;i++)
      {
        axpy_prod(sp1,m5,sp3,true);
      }
    ta2=clock();
    dt_3=ta2-ta1;
    printf("%d\n", dt_3);

    //--sparse * symmerty upper mode //...(4)
    ta1=clock();
    for(i=0;i<count;i++)
      {
        axpy_prod(sp1,sy1,sp3,true);
      }
    ta2=clock();
    dt_4=ta2-ta1;
    printf("%d\n", dt_4);

    //--sparse * symmerty lower mode //...(5)
    ta1=clock();
    for(i=0;i<count;i++)
      {
        axpy_prod(sp1,sy2,sp3,true);
      }
    ta2=clock();
    dt_5=ta2-ta1;
    printf("%d\n", dt_5);

    //--dense * dense  //...(5)
    {
      int ta1=clock();
      for(i=0;i<count;i++)
        {
          axpy_prod(m4,m5,m6,true);
        }
      int ta2=clock();
      dt_5=ta2-ta1;
      printf("%d\n", dt_5);
    }
    //--sparse * sparse column mode //...(1)
    if (0) {
      tb1=clock();
      for(i=0;i<count;i++)
        {
          axpy_prod(sp1,sp2,m6,true);
        }
      tb2=clock();
      dt_1=tb2-tb1;   
      printf("%d\n", dt_1);
    }
  }

  return 0;
}