int main(int argc, char **argv) {
    srand48(time(NULL));

    int m = atoi(argv[1]);
    int k = atoi(argv[2]);
    int n = atoi(argv[3]);
    // std::string interleaving = argv[4];
    char* outfile = argv[4];

    FILE *f = fopen(outfile,"a");

    float *A, *B, *C;
    A = (float*) malloc(m * k * sizeof(float));
    B = (float*) malloc(k * n * sizeof(float));
    C = (float*) malloc(m * n * sizeof(float));
    initialize(m, k, n, A, B, C);

    // Time multiplication
    struct timeval start, end;
    gettimeofday(&start, NULL);
    cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, m,n,k, 1, A,m, B,k, 1, C,m);
    gettimeofday(&end, NULL);
    double seconds = (end.tv_sec - start.tv_sec) + 1.0e-6 * (end.tv_usec - start.tv_usec);
    double Gflop_s = 2e-9 * m * k * n / seconds;
    fprintf(f,"MKL-SINGLE,%d,%d,%d,%s,%f\n", m, k, n, "MKL", Gflop_s);
    printf("MKL-SINGLE,%d,%d,%d,%s,%f\n", m, k, n, "MKL", Gflop_s);

    // Housekeeping
    free(A);
    free(B);
    free(C);
    fclose(f);
    return 0;
}

/******************************
compute a corr-matrix based on trial, starting row and step
input: a trial struct, starting row id, step (the row of the correlation matrix, whose column is row), the raw matrix struct array
output: the corr-matrix struct
*******************************/
CorrMatrix* CorrMatrixComputation(Trial trial, int sr, int step, RawMatrix** matrices1, RawMatrix** matrices2)
{
  int sid = trial.sid;
  int sc = trial.sc;
  int ec = trial.ec;
  int row1 = matrices1[sid]->row;
  int row2 = matrices2[sid]->row;
  int col = matrices1[sid]->col;  // the column of 1 and 2 should be the same, i.e. the number of TRs of a block
  float* mat1 = matrices1[sid]->matrix;
  float* mat2 = matrices2[sid]->matrix;
  float* buf1 = new float[row1*col]; // col is more than what really need, just in case
  float* buf2 = new float[row2*col]; // col is more than what really need, just in case
  int ml1 = getBuf(sc, ec, row1, col, mat1, buf1);  // get the normalized matrix, return the length of time points to be computed
  int ml2 = getBuf(sc, ec, row2, col, mat2, buf2);  // get the normalized matrix, return the length of time points to be computed, m1==m2
  CorrMatrix* c_matrix = new CorrMatrix();
  c_matrix->sid = sid;
  c_matrix->tlabel = trial.label;
  c_matrix->sr = sr;
  c_matrix->step = step;
  c_matrix->nVoxels = row2; //
  float* corrs = new float[step*row2];
  cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, step, row2, ml1, 1.0, buf1+sr*ml1, ml1, buf2, ml2, 0.0, corrs, row2);
  c_matrix->matrix = corrs;
  delete[] buf1;
  delete[] buf2;
  return c_matrix;
}

// Calculate distances with matrix/matrix operations (BLAS3)
void distance(int N, int D, float *data, float *result) {
	int blockSize = 512;
  int blockCount = N / blockSize;
  int remainder = N % blockSize;
  
  if (remainder) blockCount++; // Include the ragged edge
  
  dispatch_queue_t queue = dispatch_get_global_queue(DISPATCH_QUEUE_PRIORITY_DEFAULT, 0);
  
	float *diag = (float *)malloc(sizeof(float) * N);
  float *C = (float *)malloc(blockCount * blockSize * blockSize * sizeof(float));
  
  dispatch_apply(blockCount, queue, ^(size_t m) {
    int i, j;
    
    int outerDim = blockSize;
    if (m == blockCount - 1 && remainder) outerDim = remainder;
    
    
    cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, outerDim, outerDim, D,
                1, &data[m * blockSize * D], D, &data[m * blockSize * D], D, 0,
                &C[m * blockSize * blockSize], outerDim);
    
    for(i = 0; i < outerDim; i++)
      diag[m * blockSize + i] = C[m * blockSize * blockSize + i * (outerDim + 1)];
    
    for(i = 0; i < outerDim; i++)
      for(j = i + 1; j < outerDim; j++)
        result[utndidx(i + m * blockSize, j + m * blockSize)] = \
        sqrt(diag[i + m * blockSize] + diag[j + m * blockSize] - \
             2 * C[m * blockSize * blockSize + j * outerDim + i]);
  });

void col_major_fmatrix_multiply(__CLPK_integer row1_ct, __CLPK_integer col2_ct, __CLPK_integer common_ct, float* inmatrix1, float* inmatrix2, float* outmatrix) {
#ifdef NOLAPACK
    uintptr_t row1_ct_l = row1_ct;
    uintptr_t col2_ct_l = col2_ct;
    uintptr_t common_ct_l = common_ct;
    uintptr_t row_idx;
    uintptr_t col_idx;
    uintptr_t com_idx;
    float* fptr;
    float fxx;
    // not optimized
    for (col_idx = 0; col_idx < col2_ct_l; col_idx++) {
        for (row_idx = 0; row_idx < row1_ct_l; row_idx++) {
            fxx = 0;
            fptr = &(inmatrix2[col_idx * common_ct]);
            for (com_idx = 0; com_idx < common_ct_l; com_idx++) {
                fxx += (*fptr++) * inmatrix1[com_idx * row1_ct_l + row_idx];
            }
            *outmatrix++ = fxx;
        }
    }
#else
#ifdef _WIN32
    char blas_char = 'N';
    float fyy = 1;
    float fzz = 0;
    sgemm_(&blas_char, &blas_char, &row1_ct, &col2_ct, &common_ct, &fyy, inmatrix1, &row1_ct, inmatrix2, &common_ct, &fzz, outmatrix, &row1_ct);
#else
    cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, row1_ct, col2_ct, common_ct, 1.0, inmatrix1, row1_ct, inmatrix2, common_ct, 0.0, outmatrix, row1_ct);
#endif // _WIN32
#endif // NOLAPACK
}

/* C <- alpha*A*B + beta*C */
void gemm(double alpha, Mat mA, Mat mB, double beta, Mat mC) {
  const int n = MatN(mA);
  const void* const a = MatElems(mA);
  const void* const b = MatElems(mB);
  void* const c = MatElems(mC);
  const bool dev = MatDev(mA);

  switch (MatElemSize(mA)) {
  case 4:
    if (dev) {
      float alpha32 = alpha, beta32 = beta;
      cublasSgemm(g_cublasHandle, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n,
                  &alpha32, a, n, b, n, &beta32, c, n);
    } else {
      cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n,
                  n, alpha, a, n, b, n, beta, c, n);
    }
    break;

  case 8:
    if (dev) {
      cublasDgemm(g_cublasHandle, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n,
                  &alpha, a, n, b, n, &beta, c, n);
    } else {
      cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n,
                  n, alpha, a, n, b, n, beta, c, n);
    }
    break;
  }
}

void SENNA_nn_temporal_convolution(float *output, int output_frame_size,
                                   float *weights, float *biases, float *input,
                                   int input_frame_size, int n_frames,
                                   int k_w) {
#ifdef USE_BLAS
  if (k_w == 1) {
    if (biases) {
      int t;
      for (t = 0; t < n_frames; t++)
        cblas_scopy(output_frame_size, biases, 1,
                    output + t * output_frame_size, 1);
    }
    cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, output_frame_size,
                n_frames, input_frame_size, 1.0, weights, input_frame_size,
                input, input_frame_size, (biases ? 1.0 : 0.0), output,
                output_frame_size);
  } else
#endif
  {
    int t;

    for (t = 0; t < n_frames - k_w + 1; t++)
      SENNA_nn_linear(output + t * output_frame_size, output_frame_size,
                      weights, biases, input + t * input_frame_size,
                      input_frame_size * k_w);
  }
}

void ccv_gemm(ccv_matrix_t* a, ccv_matrix_t* b, double alpha, ccv_matrix_t* c, double beta, int transpose, ccv_matrix_t** d, int type)
{
	ccv_dense_matrix_t* da = ccv_get_dense_matrix(a);
	ccv_dense_matrix_t* db = ccv_get_dense_matrix(b);
	ccv_dense_matrix_t* dc = (c == 0) ? 0 : ccv_get_dense_matrix(c);

	assert(CCV_GET_DATA_TYPE(da->type) == CCV_GET_DATA_TYPE(db->type) && CCV_GET_CHANNEL(da->type) == 1 && CCV_GET_CHANNEL(db->type) == 1 && ((transpose & CCV_A_TRANSPOSE) ? da->rows : da->cols) == ((transpose & CCV_B_TRANSPOSE) ? db->cols : db->rows));

	if (dc != 0)
		assert(CCV_GET_DATA_TYPE(dc->type) == CCV_GET_DATA_TYPE(da->type) && CCV_GET_CHANNEL(dc->type) == 1 && ((transpose & CCV_A_TRANSPOSE) ? da->cols : da->rows) == dc->rows && ((transpose & CCV_B_TRANSPOSE) ? db->rows : db->cols) == dc->cols);

	ccv_declare_derived_signature_case(sig, ccv_sign_with_format(20, "ccv_gemm(%d)", transpose), ccv_sign_if(dc == 0 && da->sig != 0 && db->sig != 0, da->sig, db->sig, CCV_EOF_SIGN), ccv_sign_if(dc != 0 && da->sig != 0 && db->sig != 0 && dc->sig != 0, da->sig, db->sig, dc->sig, CCV_EOF_SIGN));
	type = CCV_GET_DATA_TYPE(da->type) | CCV_GET_CHANNEL(da->type);
	ccv_dense_matrix_t* dd = *d = ccv_dense_matrix_renew(*d, (transpose & CCV_A_TRANSPOSE) ? da->cols : da->rows, (transpose & CCV_B_TRANSPOSE) ? db->rows : db->cols, type, type, sig);
	ccv_object_return_if_cached(, dd);

	if (dd != dc && dc != 0)
		memcpy(dd->data.u8, dc->data.u8, dc->step * dc->rows);
	else if (dc == 0) // clean up dd if dc is not provided
		memset(dd->data.u8, 0, dd->step * dd->rows);

#if (defined HAVE_CBLAS || defined HAVE_ACCELERATE_FRAMEWORK)
	switch (CCV_GET_DATA_TYPE(dd->type))
	{
		case CCV_32F:
			cblas_sgemm(CblasRowMajor, (transpose & CCV_A_TRANSPOSE) ? CblasTrans : CblasNoTrans, (transpose & CCV_B_TRANSPOSE) ? CblasTrans : CblasNoTrans, dd->rows, dd->cols, (transpose & CCV_A_TRANSPOSE) ? da->rows : da->cols, alpha, da->data.f32, da->cols, db->data.f32, db->cols, beta, dd->data.f32, dd->cols);
			break;
		case CCV_64F:
			cblas_dgemm(CblasRowMajor, (transpose & CCV_A_TRANSPOSE) ? CblasTrans : CblasNoTrans, (transpose & CCV_B_TRANSPOSE) ? CblasTrans : CblasNoTrans, dd->rows, dd->cols, (transpose & CCV_A_TRANSPOSE) ? da->rows : da->cols, alpha, da->data.f64, da->cols, db->data.f64, db->cols, beta, dd->data.f64, dd->cols);
			break;
	}
#else
	assert(0 && "You need a BLAS compatible library for this function, e.g. libatlas.");
#endif
}

/***********************************************
Get the inner product of vectors from start row(sr), last rowLength-length
input: the number of subjects, the number of blocks, the start row, the number of voxels of masked matrix one that involved in the computing, the trials information, the first masked data array, the second masked data array
output: the partial similarity matrix based on the selected rows of first matrices and the whole second matrices
************************************************/
float* GetPartialInnerSimMatrixWithMasks(int nSubs, int nTrials, int sr, int rowLength, Trial* trials, RawMatrix** masked_matrices1, RawMatrix** masked_matrices2) // compute the correlation between masked matrices
{
  int i;
  int row1 = masked_matrices1[0]->row;
  int row2 = masked_matrices2[0]->row;  //rows should be the same across subjects since we are using the same mask file to filter out voxels
  float* values= new float[nTrials*rowLength*row2];
  float* simMatrix = new float[nTrials*nTrials];
  memset((void*)simMatrix, 0, nTrials*nTrials*sizeof(float));
  for (i=0; i<nTrials; i++)
  {
    int sc = trials[i].sc;
    int ec = trials[i].ec;
    int sid = trials[i].sid;
    int col = masked_matrices1[sid]->col; // the column of 1 and 2 should be the same, i.e. the number of TRs of a block; columns may be different, since different subjects have different TRs
    float* mat1 = masked_matrices1[sid]->matrix;
    float* mat2 = masked_matrices2[sid]->matrix;
    float* buf1 = new float[row1*col]; // col is more than what really need, just in case
    float* buf2 = new float[row2*col]; // col is more than what really need, just in case
    int ml1 = getBuf(sc, ec, row1, col, mat1, buf1);  // get the normalized matrix, return the length of time points to be computed
    int ml2 = getBuf(sc, ec, row2, col, mat2, buf2);  // get the normalized matrix, return the length of time points to be computed, m1==m2
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, rowLength, row2, ml1, 1.0, buf1+sr*ml1, ml1, buf2, ml2, 0.0, values+i*rowLength*row2, row2);
    delete[] buf1;
    delete[] buf2;
  }
  NormalizeCorrValues(values, nTrials, rowLength, row2, nSubs);
  GetDotProductUsingMatMul(simMatrix, values, nTrials, rowLength, row2);
  delete[] values;
  return simMatrix;
}

// row here is nTops, get the inner product of vectors from start row(sr), last rowLength-length
float* GetPartialInnerSimMatrix(int row, int col, int nSubs, int nTrials, int sr, int rowLength, Trial* trials, RawMatrix** r_matrices) // only compute the correlation among the selected voxels
{
  int i;
  float* values = new float[nTrials*rowLength*row];
  float* simMatrix = new float[nTrials*nTrials];
  for (i=0; i<nTrials*nTrials; i++) simMatrix[i] = 0.0;
  for (i=0; i<nTrials; i++)
  {
    int sc = trials[i].sc;
    int ec = trials[i].ec;
    int sid = trials[i].sid;
    float* mat = r_matrices[sid]->matrix;
    //if (i==0 && sr==0) cout<<mat[1000*col]<<" "<<mat[1000*col+1]<<" "<<mat[1000*col+2]<<" "<<mat[1000*col+3]<<endl;
    //else if (i==0 && sr!=0) cout<<mat[0]<<" "<<mat[1]<<" "<<mat[2]<<" "<<mat[3]<<endl;
    float* buf = new float[row*col]; // col is more than what really need, just in case
    int ml = getBuf(sc, ec, row, col, mat, buf);  // get the normalized matrix, return the length of time points to be computed
    //cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, step, row, ml, 1.0, buf+sr*ml, ml, buf, ml, 0.0, corrs, row);
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, rowLength, row, ml, 1.0, buf+sr*ml, ml, buf, ml, 0.0, values+i*rowLength*row, row);
    delete[] buf;
  }
  NormalizeCorrValues(values, nTrials, rowLength, row, nSubs);
  GetDotProductUsingMatMul(simMatrix, values, nTrials, rowLength, row);
  delete[] values;
  return simMatrix;
}

int
main (int argc, char **argv) {
	// Log messeages into stderr
	FLAGS_logtostderr = 1;
	google::InitGoogleLogging (argv[0]);

	LOG(INFO) << "Begin BLAS Demo";
	std::cout << "  A = [ 0 0 1 ; 0 1 0 ; 1 0 0 ]" << std::endl
		<< "  B = [ 1;2;3 ]" << std::endl;
	int M=3, K=3, N=1;
	float A[M*K] = { 0,0,1, 0,1,0, 1,0,0 };
	float B[K*N] = { 1,2,3 };
	float C[M*N];

	// Do matrix multiplication
	LOG(INFO) << "Calculating  C := 1.0 A * B + 0.0 C ";
	cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, M, N, K,
		1.0, A, K, B, N, 0.0, C, N);

	// Dump matrix	
	LOG(INFO) << "Dumping matrix C";
	std::cout << "C = [" << std::endl;
	for (int i = 1; i <= M; i++) {
		std::cout << "\t";
		for (int j = 1; j <= N; j++) {
			std::cout << C[i*j-1] << " ";
		}
		std::cout << std::endl;
	}
	std::cout << "];" << std::endl;

	LOG(INFO) << "Demo done.";
	return 0;
}

inline void gemm( const Order order, const TransA transa, const TransB transb,
        const int m, const int n, const int k, const float alpha,
        const float* a, const int lda, const float* b, const int ldb,
        const float beta, float* c, const int ldc ) {
    cblas_sgemm( cblas_option< Order >::value, cblas_option< TransA >::value,
            cblas_option< TransB >::value, m, n, k, alpha, a, lda, b, ldb,
            beta, c, ldc );
}

void gemm(bool transa, bool transb, int m, int n, int k, float alpha, const float* A, int lda,
    const float* B, int ldb, float beta, float* C, int ldc)
{
  const CBLAS_TRANSPOSE ctransa = transa ? CblasTrans : CblasNoTrans;
  const CBLAS_TRANSPOSE ctransb = transb ? CblasTrans : CblasNoTrans;

  cblas_sgemm(CblasColMajor, ctransa, ctransb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
}

void caffe_cpu_gemm<float>(const CBLAS_TRANSPOSE TransA,
                           const CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
                           const float alpha, const float* A, const float* B, const float beta,
                           float* C) {
    int lda = (TransA == CblasNoTrans) ? K : M;
    int ldb = (TransB == CblasNoTrans) ? N : K;
    cblas_sgemm(CblasRowMajor, TransA, TransB, M, N, K, alpha, A, lda, B,
                ldb, beta, C, N);
}

inline void STARPU_SGEMM(char *transa, char *transb, int M, int N, int K, 
			float alpha, const float *A, int lda, const float *B, int ldb, 
			float beta, float *C, int ldc)
{
	enum CBLAS_TRANSPOSE ta = (toupper(transa[0]) == 'N')?CblasNoTrans:CblasTrans;
	enum CBLAS_TRANSPOSE tb = (toupper(transb[0]) == 'N')?CblasNoTrans:CblasTrans;

	cblas_sgemm(CblasColMajor, ta, tb,
			M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);				
}

void caffe_cpu_gemm<float>(const CBLAS_TRANSPOSE TransA,
    const CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
    const float alpha, const float* A, const float* B, const float beta,
    float* C) {
  int lda = (TransA == CblasNoTrans) ? K : M;
  int ldb = (TransB == CblasNoTrans) ? N : K;
#if 0
  LOG(INFO)<<"\t\t----> XEON: M="<< M <<" N="<< N <<" K="<< K;
#endif
  cblas_sgemm(CblasRowMajor, TransA, TransB, M, N, K, alpha, A, lda, B,
      ldb, beta, C, N);
}

 inline
 void gemm (CBLAS_ORDER const Order,
            CBLAS_TRANSPOSE const TransA, CBLAS_TRANSPOSE const TransB,
            int const M, int const N, int const K,
            float const alpha, float const* A, int const lda,
            float const* B, int const ldb,
            float const beta, float* C, int const ldc)
 {
   cblas_sgemm (Order, TransA, TransB, M, N, K,
                alpha, A, lda,
                B, ldb,
                beta, C, ldc);
 }

void DenseTransitionMatrix<float>::mult(const DenseTransitionMatrix<float>& A,
		const DenseTransitionMatrix<float>& B) {

	assert(A.n == B.n && A.n == n);
	assert(A.ld == B.ld && A.ld == ld);

	const float alpha_d = 1.0;
	const float beta_d = 0.0;

	// use cblas
	cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, A.n, A.n, A.n,
			alpha_d, B.data, B.ld, A.data, A.ld, beta_d, data, ld);
}

Matrix& Matrix::operator*=(const Matrix& m) throw()
{
	assert(col == m.row);
	Matrix temp(*this);
	if(col != m.col)
	{
		delete [] num;
		col = m.col;
		num = new float[row*col];
	}
	cblas_sgemm(R, N, N, row, col, m.row, 1.0f, 
	temp.num, m.row, m.num, m.col, 0.0f, num, col);
	return *this;
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_sgemm 
(JNIEnv * env, jobject calling_obj, jint order, jint transA, jint transB, jint M, jint N, jint K, 
jfloat alpha, jfloatArray jA, jint lda, jfloatArray jB, jint ldb, jfloat beta, jfloatArray jC, jint ldc){
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
	jfloat * B = (*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE);
	jfloat * C = (*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE);

	cblas_sgemm((CBLAS_ORDER)order, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K, 
                    alpha, A, lda, B, ldb, beta, C, ldc);

	(*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

void inner_multiply(int M, int K, int m, int k, int n, float *A, float *B, float *C, int depth, int CM) {
  if (depth >= MAX_DEPTH) {
    cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, m,n,k, 1, A,M, B,K, 0, C, CM);
    return;
  }

  int next_depth = depth + 1;

  int max = k;
  if (m > max) {
    max = m;
  }
  if (n > max) {
    max = n;
  }

  if (max == n) {
    n = n/2;
    float *B1 = B;
    float *B2 = B + n*K;
    cilk_spawn inner_multiply(M, K, m, k, n, A, B1, C, next_depth, CM);
    inner_multiply(M, K, m, k, n, A, B2, C + n*CM, next_depth, CM);
    cilk_sync;

  } else if (max == m) {
    m = m/2;
    float *A1 = A;
    float *A2 = A + m;
    cilk_spawn inner_multiply(M, K, m, k, n, A1, B, C, next_depth, CM);
    inner_multiply(M, K, m, k, n, A2, B, C + m, next_depth,CM);
    cilk_sync;

  } else {
    k = k/2;
    float *A1 = A;
    float *A2 = A + k*M;
    float *B1 = B;
    float *B2 = B + k;
    float *Q1 = (float*) malloc(m * n * sizeof(float));
    cilk_spawn inner_multiply(M, K, m, k, n, A1, B1, Q1, next_depth, m);
    inner_multiply(M, K, m, k, n, A2, B2, C, next_depth, CM);
    cilk_sync;
    int x;
    for (x = 0; x < n; x++) {
      cblas_saxpy(m, 1, Q1 + m*x, 1, C + CM*x, 1);
    }
    free(Q1);
  }
}