Initial modification

analysis.R

+library(ggplot2)
+library(reshape2)
+custom_scal = t(read.table("custom_scal.txt", sep=",", skipNul = T)[1:10])
+eigen_scal = t(read.table("eigen_scal.txt", sep=",", skipNul = T)[1:10])
+matrix_sizes = 3:8
+
+custom_scal = data.frame('Custom', custom_scal)
+eigen_scal = data.frame('Eigen', eigen_scal)
+colnames(custom_scal) = c('program', matrix_sizes)
+colnames(eigen_scal) = c('program', matrix_sizes)
+
+measures = rbind(melt(custom_scal), melt(eigen_scal))
+
+linear_scal = mean(custom_scal$`3`)/9*(3:8)**2
+
+ggplot(measures, aes(variable, value, group=program, color=program)) +
+  geom_smooth(se = FALSE) +
+  geom_point() +
+  xlab("Matrix Size") +
+  ylab("Time (s)") +
+  ggtitle("Custom implementation vs Eigen time comparison (2e8 matrix multiplications)") +
+  theme_minimal()
+
+ggplot(measures, aes(as.numeric(levels(variable))[variable]**2, value, group=program, color=program)) +
+  geom_smooth(se = FALSE) +
+  geom_point() +
+  xlab("Number of elements") +
+  ylab("Time (s)") +
+  ggtitle("Custom implementation vs Eigen time comparison (2e8 matrix multiplications)") +
+  theme_minimal()
+
+speedup = colMeans(eigen_scal[2:7])/colMeans(custom_scal[2:7])
+
+ggplot(stack(speedup), aes(ind, values)) +
+  geom_bar(stat="identity") +
+  xlab("Matrix Size") +
+  ylab("Speedup (Eigen time / Custom time)") +
+  ggtitle("Custom implementation vs Eigen speedup") +
+  theme_minimal() +
+  geom_hline(yintercept=1, color="red") +
+  geom_text(aes(ind, values + 0.2, label=round(values, digits=1)))
+
+
+ggplot(measures, aes(as.numeric(levels(variable))[variable]**2, value/(as.numeric(levels(variable))[variable]**2), group=program, color=program)) +
+  geom_smooth(se = FALSE) +
+  geom_point() +
+  xlab("Number of elements") +
+  ylab("Time per element (s)") +
+  ggtitle("Custom implementation vs Eigen time per element comparison") +
+  theme_minimal()
\ No newline at end of file

custom_scal.txt

+3.86, 3.86, 3.85, 3.85, 3.85, 3.87, 3.87, 3.86, 3.86, 3.85, 
+4.50, 4.50, 4.46, 4.52, 4.52, 4.50, 4.49, 4.46, 4.50, 4.50, 
+5.47, 5.40, 5.37, 5.65, 5.41, 5.40, 5.38, 5.45, 5.66, 5.63, 
+7.74, 7.94, 7.72, 7.94, 7.92, 7.94, 7.74, 7.92, 7.95, 7.72, 
+10.72, 10.68, 10.88, 10.95, 11.20, 10.68, 10.99, 10.68, 10.70, 10.67, 
+14.21, 14.13, 14.42, 14.49, 14.54, 14.21, 14.44, 14.65, 14.54, 14.22, 

eigen_scal.txt

+7.04, 7.04, 7.02, 7.01, 7.03, 7.08, 7.06, 7.01, 7.08, 7.01, 
+6.53, 6.54, 6.47, 6.48, 6.50, 6.50, 6.52, 6.50, 6.47, 6.50, 
+9.08, 9.14, 9.10, 9.08, 8.99, 9.22, 9.13, 9.15, 9.14, 9.16, 
+20.24, 20.20, 20.09, 19.56, 20.11, 19.93, 19.98, 19.99, 19.99, 19.98, 
+54.16, 54.40, 54.20, 54.49, 54.42, 54.65, 54.90, 54.67, 55.21, 54.30, 
+81.90, 81.77, 82.11, 81.57, 81.68, 85.79, 82.81, 82.43, 81.81, 84.20, 

matrix_product_custom.cpp

+#include<boost/simd/pack.hpp>
+#include<boost/simd/constant/zero.hpp>
+#include<boost/simd/function/aligned_store.hpp>
+#include<iostream>
+#include<cstring>
+
+
+#ifndef SIZE
+#define SIZE 4
+#endif
+
+namespace bs = boost::simd;
+
+template <typename T, int MatSize, int VecSize>
+class PaddedMatrix;
+template <typename T, int MatSize, int VecSize>
+void matrix_mul(PaddedMatrix<T, MatSize, VecSize> &a,
+                PaddedMatrix<T, MatSize, VecSize> &b,
+                PaddedMatrix<T, MatSize, VecSize> &c);
+
+template <typename T, int MatSize, int VecSize>
+std::ostream & operator<<(std::ostream &os, PaddedMatrix<T, MatSize, VecSize> &mat) {
+    for (int i=0; i<MatSize; i++) {
+        for (int j=0; j<MatSize; j++) {
+            os << mat.array[i*VecSize + j] << ", ";
+        }
+        os << std::endl;
+    }
+    return os;
+}
+
+template <typename T, int MatSize, int VecSize>
+class PaddedMatrix {
+    static_assert(MatSize <= VecSize, "Matrix size is greater than vector size");
+    using matrix_t = PaddedMatrix<T, MatSize, VecSize>;
+
+    public:
+    matrix_t operator=(float const a[]) {
+        if (MatSize == VecSize) {
+            std::memcpy(this->array, a, sizeof(T)*MatSize*VecSize);
+        }
+        else {
+            for (int i=0; i<MatSize; i++) {
+                std::memcpy(&this->array[i*VecSize], &a[i*MatSize], sizeof(T)*MatSize);
+            }
+        }
+        return *this;
+    }
+    matrix_t operator=(PaddedMatrix<T, MatSize, VecSize> &other) {
+        std::memcpy(this->array, other.array, sizeof(T)*MatSize*VecSize);
+        return *this;
+    }
+
+
+    friend void matrix_mul<T, MatSize, VecSize>(matrix_t &a, matrix_t &b, matrix_t &c);
+    friend std::ostream & operator<<<T, MatSize, VecSize>(std::ostream &os, matrix_t &mat);
+    private:
+    alignas(sizeof(T)*VecSize) float array[MatSize*VecSize] = {0};
+
+};
+
+template <typename T, int MatSize, int VecSize>
+static inline void matrix_mul(PaddedMatrix<T, MatSize, VecSize> &a,
+                PaddedMatrix<T, MatSize, VecSize> &b,
+                PaddedMatrix<T, MatSize, VecSize> &c)
+{
+    using pack_t = bs::pack<T, VecSize>;
+    for (int i=0; i<MatSize; i++) {
+        pack_t res = bs::Zero<pack_t>();
+        int mult = VecSize*i;
+        for (int j=0; j<MatSize; j++) {
+            pack_t factor {a.array[mult + j]};
+            pack_t row(&b.array[VecSize*j]);
+            res += factor * row;
+        }
+        bs::aligned_store(res, &c.array[mult]);
+    }
+}
+
+int main(int argc, char **argv) {
+    alignas(32) float a[SIZE*SIZE];
+    alignas(32) float b[SIZE*SIZE];
+
+    for (int i=0; i<SIZE; i++) {
+        for (int j=0; j<SIZE; j++) {
+            a[i*SIZE+j] = i+j;
+        }
+    }
+
+    for (int i=0; i<SIZE; i++) {
+        for (int j=0; j<SIZE; j++) {
+            float val;
+            if (i == 0 && j == 1) val = -1;
+            else if (i == 1 && j == 0) val = 1;
+            else if (i > 1 && i == j && i % 2) val = -1;
+            else if (i > 1 && i == j && !(i % 2)) val = 1;
+            else val = 0;
+            b[i*SIZE+j] = val;
+        }
+    }
+
+    PaddedMatrix<float, SIZE, 8> pa, pb, pf;
+    pa = a;
+    pb = b;
+
+//    std::cout << pa << std::endl;
+//    std::cout << pb << std::endl;
+
+    for (unsigned int i=0; i<200000000; i++) {
+        matrix_mul<float, SIZE, 8>(pa, pb, pf);
+        matrix_mul<float, SIZE, 8>(pf, pb, pa);
+    }
+
+    std::cout << pa << std::endl;
+    return 0;
+}

matrix_product_eigen.cpp

+#include<eigen3/Eigen/Dense>
+#include<iostream>
+
+
+#ifndef SIZE
+#define SIZE 4
+#endif
+
+int main(int argc, char **argv) {
+    Eigen::Matrix<float, SIZE, SIZE> pa, pb, pf;
+
+    for (int i=0; i<SIZE; i++) {
+        for (int j=0; j<SIZE; j++) {
+            pa(i,j) = i+j;
+        }
+    }
+
+    for (int i=0; i<SIZE; i++) {
+        for (int j=0; j<SIZE; j++) {
+            float val;
+            if (i == 0 && j == 1) val = -1;
+            else if (i == 1 && j == 0) val = 1;
+            else if (i > 1 && i == j && i % 2) val = -1;
+            else if (i > 1 && i == j && !(i % 2)) val = 1;
+            else val = 0;
+            pb(i,j) = val;
+        }
+    }
+
+//    std::cout << pa << std::endl;
+//    std::cout << pb << std::endl;
+
+    for (unsigned int i=0; i<200000000; i++) {
+        pf = pb * pa;
+        pa = pb * pf;
+    }
+
+    std::cout << pa << std::endl;
+    return 0;
+}

measure_perf.sh

+#!/bin/bash
+export REPEAT=10
+export PROG="a.out"
+export FILE="scal.txt"
+export TIME=/usr/bin/time
+export BOOST_SIMD_ROOT=/home/ping/These/boost.simd/include
+
+rm -f custom_${FILE}
+rm -f eigen_${FILE}
+
+for size in $(seq 3 8); do
+    g++ matrix_product_custom.cpp -I ${BOOST_SIMD_ROOT} -DSIZE=$size -O3 -march=native
+    for i in $(seq 1 ${REPEAT}); do
+        /usr/bin/time -f "%U, " ./a.out 2> >(tr -d '\n' | tee -a custom_${FILE}) && echo
+    done
+    echo >> custom_${FILE}
+
+    g++ matrix_product_eigen.cpp -DSIZE=$size -O3 -march=native
+    for i in $(seq 1 ${REPEAT}); do
+        /usr/bin/time -f "%U, " ./a.out 2> >(tr -d '\n' | tee -a eigen_${FILE}) && echo
+    done
+    echo >> eigen_${FILE}
+done