#include "../mat_checks.h"
#include <itpp/itbase.h>
using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

//#define NDEBUG

mat matmul ( mat &A, mat &B ) {
    mat C ( A.rows(), B.cols() );
    int i, j, k;
    double sum;

    for ( i = 0; i < A.rows(); i++ ) {
        for ( j = 0; j < A.cols(); j++ ) {
            sum = 0.0;
            for ( k = 0; k < A.cols(); k++ ) {
                sum += A._elem ( i, k ) * B._elem ( k, j );
            }
            C ( i, j ) = sum;
        }
    }
    return C;
}

void matmul2 ( int n,  double *A, double *B, double *C ) {
    int i, j, k;
    double sum;

    for ( i = 0; i < n; i++ ) {
        for ( j = 0; j < n; j++ ) {
            sum = 0.0;
            for ( k = 0; k < n; k++ ) {
                sum += A [ i*n+k ] * B [ k*n+j ];
            }
            C[ i*n+j] = sum;
        }
    }
//	return C;
}

TEST ( blas_stress ) {
    Real_Timer tt;
    vec exec_times ( 4 );
    vec exec_times_b ( 4 );
    vec exec_times_c ( 4 );

    mat A;
    mat B;
    mat C;

    vec vn = "5 50 200 500";
    int n;

    for ( int i = 0; i < vn.length(); i++ ) {
        n = (int) vn ( i );
        A = randu ( n, n );
        B = randu ( n, n );

        tt.tic();
        for ( int ii = 0; ii < 10; ii++ ) {
            C = matmul ( A, B );
        }
        exec_times ( i ) = tt.toc();

        tt.tic();
        for ( int ii = 0; ii < 10; ii++ ) {
            C = A * B;
        }
        exec_times_b ( i ) = tt.toc();

        C = zeros ( n, n );
        tt.tic();
        for ( int ii = 0; ii < 10; ii++ ) {
            matmul2 ( n, A._data(), B._data(), C._data() );
        }
        exec_times_c ( i ) = tt.toc();
    }
    cout << exec_times << endl;
    cout << exec_times_b << endl;
    cout << exec_times_c << endl;

    it_file itf ( "blas_test.it" );
    itf << Name ( "exec_times" ) << exec_times;
    itf << Name ( "exec_times_b" ) << exec_times_b;
    itf << Name ( "exec_times_c" ) << exec_times_c;
}
