From bd3662a047815902544277ae6bef9f29dea814a2 Mon Sep 17 00:00:00 2001
From: "Ma, Jiangnan" <jiangnan.ma@intel.com>
Date: Mon, 9 Oct 2023 05:34:10 +0800
Subject: [PATCH 1/2] add matmul version1

---
 high_performance_operators/gemm/version1.cpp | 82 ++++++++++++++++++++
 1 file changed, 82 insertions(+)
 create mode 100644 high_performance_operators/gemm/version1.cpp

diff --git a/high_performance_operators/gemm/version1.cpp b/high_performance_operators/gemm/version1.cpp
new file mode 100644
index 0000000..e650d3d
--- /dev/null
+++ b/high_performance_operators/gemm/version1.cpp
@@ -0,0 +1,82 @@
+#include <iostream>
+#include <vector>
+#include <cassert>
+#include <sycl/sycl.hpp>
+
+// Transposition
+
+// one should always start with memory before optimizing arithmetic, as it is much more likely to be the bottleneck.
+
+template<class type>
+struct matrix_2d {
+    matrix_2d(std::vector<std::vector<type>> &data) {
+        this->data = data;
+        h = data.size();
+        w = data[0].size();
+    }
+
+    matrix_2d(int high, int width, type fill = 0) {
+        data.resize(high, std::vector<type>(width, fill));
+        h = high;
+        w = width;
+    }
+
+    std::vector<type> & operator[](int i) {
+        return data[i];
+    }
+
+    friend std::ostream & operator<<(std::ostream & o, matrix_2d<type> & m) {
+        for (int i = 0; i < m.h; ++i) {
+            for (int j = 0; j < m.w; ++j) o << m[i][j] << " ";
+            o << std::endl;
+        }
+        return o;
+    }
+
+    int h;
+    int w;
+    std::vector<std::vector<type>> data;
+};
+
+template<class type>
+void transposition(matrix_2d<type> & old, matrix_2d<type> & trans) {
+    assert(old.h == trans.w && old.w == trans.h);
+
+    for (int i = 0; i < old.h; ++i) {
+        for (int j = 0; j < old.w; ++j) {
+            trans[j][i] = old[i][j];
+        }
+    }
+    return;
+}
+
+template<class type>
+matrix_2d<type> & gemm2d_version1(matrix_2d<type> & a, matrix_2d<type> & b, matrix_2d<type> & result) {
+    assert(a.w == b.h);
+    assert(result.h == a.h && result.w == b.w);
+
+    auto trans_b = matrix_2d<type>(b.w, b.h);
+    transposition<type>(b, trans_b);
+
+    for (int i = 0; i < result.h; ++i) {
+        for (int j = 0; j < result.w; ++j) {
+            type tmp = 0;
+            for (int k = 0; k < a.w; ++k) {
+                tmp += a[i][k] * trans_b[j][k];
+            }
+            result[i][j] = tmp;
+        }
+    }
+    return result;
+};
+
+int main() {
+    auto a = matrix_2d<float>(4, 3, 2);
+    auto b = matrix_2d<float>(3, 5, 1.5);
+    auto c = matrix_2d<float>(4, 5);
+
+    gemm2d_version1<float>(a, b, c);
+    std::cout << a << b << c << std::endl;
+
+    return 0;
+}

From 4468380ddac4ec618124062d780b37d06b37eefb Mon Sep 17 00:00:00 2001
From: "Ma, Jiangnan" <jiangnan.ma@intel.com>
Date: Fri, 27 Oct 2023 06:58:39 +0800
Subject: [PATCH 2/2] improve matrix

---
 high_performance_operators/gemm/version1.cpp | 55 ++++++++++++++------
 1 file changed, 40 insertions(+), 15 deletions(-)

diff --git a/high_performance_operators/gemm/version1.cpp b/high_performance_operators/gemm/version1.cpp
index e650d3d..2b9a068 100644
--- a/high_performance_operators/gemm/version1.cpp
+++ b/high_performance_operators/gemm/version1.cpp
@@ -9,33 +9,58 @@
 
 template<class type>
 struct matrix_2d {
-    matrix_2d(std::vector<std::vector<type>> &data) {
-        this->data = data;
-        h = data.size();
-        w = data[0].size();
+
+public:
+    matrix_2d(std::vector<std::vector<type>> data) {
+      h = data.size();
+      w = data[0].size();
+      stride[0] = w;
+      stride[1] = 1;
+      
+      this->data = new type[h * w];
+      for (int i = 0; i < h; ++i) {
+        for (int j = 0; j < w; ++j) {
+          (this->data)[i * stride[0] + j * stride[1]] = data[i][j];  
+        }
+      }
     }
 
     matrix_2d(int high, int width, type fill = 0) {
-        data.resize(high, std::vector<type>(width, fill));
-        h = high;
-        w = width;
+      h = high;
+      w = width;
+      stride[0] = w;
+      stride[1] = 1;
+
+      this->data = new type[h * w];
+      std::fill(this->data, this->data + h * w, fill);
+    }
+
+    ~matrix_2d() {
+      delete data;
     }
 
-    std::vector<type> & operator[](int i) {
-        return data[i];
+    type * data_ptr() {
+      return data;
+    }
+
+    type * operator[](int i) {
+      return (this->data) + i * stride[0];
     }
 
     friend std::ostream & operator<<(std::ostream & o, matrix_2d<type> & m) {
-        for (int i = 0; i < m.h; ++i) {
-            for (int j = 0; j < m.w; ++j) o << m[i][j] << " ";
-            o << std::endl;
-        }
-        return o;
+      for (int i = 0; i < m.h; ++i) {
+          for (int j = 0; j < m.w; ++j) o << m[i][j] << " ";
+          o << std::endl;
+      }
+      return o;
     }
 
     int h;
     int w;
-    std::vector<std::vector<type>> data;
+
+private:
+    int stride[2];
+    type * data;
 };
 
 template<class type>