#include "stdio.h"
#include "stdlib.h"
#include "f77blas.h"
#include "sys/time.h"
#include "time.h"
#include <iostream>
#include <string>
#include <sstream>

using namespace std;

extern void dgemm_(char*, char*, int*, int*,int*, double*, double*, int*, double*, int*, double*, double*, int*);

string expect = R"(16.801 18.002 18.003 16.801 15.602 22.803 )";

void setup() {
  // put your setup code here, to run once:
  // Loop counter
  int i = 0;
  // Number of rows in matrix A, number of rows in matrix C
  int m = 2;
  // Number of columns in matrix B, number of columns in matrix C
  int n = 3;
  // Number of columns in matrix A, number of rows in matrix B
  int k = 4;
  // Number of elements in matrix A
  int sizeofa = m * k;
  // Number of elements in matrix B
  int sizeofb = k * n;
  // Number of elements in matrix C
  int sizeofc = m * n;
  // Flag for whether matrix A is transposed, 'N' means no transpose
  char ta = 'N';
  // Flag for whether matrix B is transposed, 'N' means no transpose
  char tb = 'N';
  // Scalar factor alpha in matrix multiplication
  double alpha = 1.2;
  // Scalar factor beta in matrix multiplication
  double beta = 0.001;

  // Dynamically allocate memory for matrix A
  double* A = (double*)malloc(sizeof(double) * sizeofa);
  // Dynamically allocate memory for matrix B
  double* B = (double*)malloc(sizeof(double) * sizeofb);
  // Dynamically allocate memory for matrix C
  double* C = (double*)malloc(sizeof(double) * sizeofc);

  // Initialize random seed
  srand((unsigned)time(NULL));

  // Initialize elements of matrix A
  for (i=0; i<sizeofa; i++)
    // Elements of matrix A initialized as i % 3 + 1
    A[i] = i%3+1;//(rand()%100)/10.0;

  // Initialize elements of matrix B
  for (i=0; i<sizeofb; i++)
    // Elements of matrix B initialized as i % 3 + 1
    B[i] = i%3+1;//(rand()%100)/10.0;

  // Initialize elements of matrix C
  for (i=0; i<sizeofc; i++)
    // Elements of matrix C initialized as i % 3 + 1
    C[i] = i%3+1;//(rand()%100)/10.0;

  // Print matrix multiplication parameter info
  printf("m=%d,n=%d,k=%d,alpha=%lf,beta=%lf,sizeofc=%d\n",m,n,k,alpha,beta,sizeofc);
  // Call dgemm_ function for matrix multiplication: C = alpha * A * B + beta * C
  dgemm_(&ta, &tb, &m, &n, &k, &alpha, A, &m, B, &k, &beta, C, &m);

  // Print info of matrix A
  printf("This is matrix A\n\n");
  for(i=0; i < sizeofa; i++)
    // Print each element of matrix A
    printf("%lf ", A[i]);
  printf("\n");

  // Print info of matrix B
  printf("This is matrix B\n\n");
  for(i=0; i < sizeofb; i++)
    // Print each element of matrix B
    printf("%lf ", B[i]);
  printf("\n");

  // Create a stringstream object to store the output of matrix C
  stringstream ss;
  // Save current standard output buffer
  streambuf   *buffer = cout.rdbuf();
  // Redirect standard output to stringstream
  cout.rdbuf(ss.rdbuf());

  // Traverse each element of matrix C
  for(int i=0;i<sizeofc;i++)
    // Output element i of matrix C followed by a space
    {cout<<C[i]<<" ";}

  // Restore standard output buffer
  cout.rdbuf(buffer);
  // Convert stringstream content to string
  string s(ss.str());
  // Print separator line
  cout << "*********************************************************" << endl;
  // Print result prompt
  cout << "This is the result:" << endl;
  // Print computed result of matrix C
  cout << s << endl;
  // Print separator line
  cout << "*********************************************************" << endl;
  // Print reference prompt
  cout << "This is the reference:" << endl;
  // Print expected result
  cout << expect << endl;

  // Compare actual result with expected result
  if (expect == s)
    // If results match, print test passed
    cout << "{Test PASS}." << endl;
  else
    // If results differ, print test failed
    cout << "{Test FAIL}." << endl;

  // Free memory for matrix A
  free(A);
  // Free memory for matrix B
  free(B);
  // Free memory for matrix C
  free(C);

}

void loop() {
  // put your main code here, to run repeatedly:

}
