#include <iostream>
#include <numeric>
using namespace std;
#include <cblas.h>
#include <xmmintrin.h>

/*
extern "C" void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
                 const int K, const double alpha, const double *A,
                 const int lda, const double *B, const int ldb,
                 const double beta, double *C, const int ldc);
 */
int DIM = 1000;
double *A;
double *B;
double *C;
                 
int main(int argc, char *argv[]) {
    char TRANS = 'N';
    double alpha = 1.0;
    double beta = 0.0;

	srand(19702013);
	
	if (argc > 1) {
		DIM = atoi(argv[1]);
	}
	
	int size = DIM*DIM;
	
	// allocate
	A = static_cast<double *>(_mm_malloc(size * sizeof(double), 16));
	B = static_cast<double *>(_mm_malloc(size * sizeof(double), 16));
	C = static_cast<double *>(_mm_malloc(size * sizeof(double), 16));
	
	// fill
	for (int i=0; i<size; ++i) {
		A[i] = (rand() % 131) * 0.01;
		B[i] = (rand() % 131) * 0.01;
	}
	
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, DIM, DIM, DIM, 
    	alpha, A, DIM, B, DIM, beta, C, DIM);

	// check result
	double sum = accumulate(&C[0], &C[size], 0.0);
	
	cout << "sum = " << sum << endl;
	
	_mm_free(A);
	_mm_free(B);
	_mm_free(C);
    return 0;
}