From d6162af210d9d5648d33bf0fda40f773ac200df5 Mon Sep 17 00:00:00 2001
From: AlexeyAB <alexeyab84@gmail.com>
Date: Wed, 08 Aug 2018 23:31:36 +0000
Subject: [PATCH] Optimized on CPU: gemm_bin, im2col, activation, transpose

---
 src/network.c             |    3 
 src/gemm.h                |   13 
 src/convolutional_layer.c |   78 ++++++
 src/gemm.c                |  466 ++++++++++++++++++++++++++++++++++++----------
 src/activations.c         |   14 +
 src/convolutional_layer.h |    2 
 6 files changed, 450 insertions(+), 126 deletions(-)

diff --git a/src/activations.c b/src/activations.c
index 0cbb2f5..eab4e23 100644
--- a/src/activations.c
+++ b/src/activations.c
@@ -95,8 +95,16 @@
 void activate_array(float *x, const int n, const ACTIVATION a)
 {
     int i;
-    for(i = 0; i < n; ++i){
-        x[i] = activate(x[i], a);
+    if (a == LINEAR) {}
+    else if (a == LEAKY) {
+        for (i = 0; i < n; ++i) {
+            x[i] = leaky_activate(x[i]);
+        }
+    }
+    else {
+        for (i = 0; i < n; ++i) {
+            x[i] = activate(x[i], a);
+        }
     }
 }
 
@@ -139,5 +147,5 @@
     for(i = 0; i < n; ++i){
         delta[i] *= gradient(x[i], a);
     }
-} 
+}
 
diff --git a/src/convolutional_layer.c b/src/convolutional_layer.c
index 0bde97a..a820588 100644
--- a/src/convolutional_layer.c
+++ b/src/convolutional_layer.c
@@ -593,15 +593,15 @@
     }
 }
 
-void binary_transpose_align_weights(convolutional_layer *l, size_t ldb_align)
+void binary_align_weights(convolutional_layer *l, size_t lda_align)
 {
     int m = l->n;
     int k = l->size*l->size*l->c;
-    size_t new_ldb = k + (ldb_align - k%ldb_align); // (k / 8 + 1) * 8;
+    size_t new_lda = k + (lda_align - k%lda_align); // (k / 8 + 1) * 8;
 
     binarize_weights(l->weights, m, k, l->binary_weights);
 
-    size_t align_weights_size = new_ldb * m;
+    size_t align_weights_size = new_lda * m;
     size_t align_bit_weights_size = align_weights_size / 8;// +1;
     float *align_weights = calloc(align_weights_size, sizeof(float));
     l->align_bit_weights = calloc(align_bit_weights_size, sizeof(char));
@@ -610,7 +610,7 @@
     // align A without transpose
     for (i = 0; i < m; ++i) {
         for (j = 0; j < k; ++j) {
-            align_weights[i*new_ldb + j] = l->binary_weights[i*k + j];
+            align_weights[i*new_lda + j] = l->binary_weights[i*k + j];
         }
     }
     float_to_bit(align_weights, l->align_bit_weights, align_weights_size);
@@ -622,6 +622,56 @@
 }
 
 
+size_t binary_transpose_align_input(int k, int n, float *b, char **t_bit_input, size_t ldb_align)
+{
+    size_t new_ldb = k + (ldb_align - k%ldb_align); // (k / 8 + 1) * 8;
+    size_t t_intput_size = new_ldb * n;
+    size_t t_bit_input_size = t_intput_size / 8;// +1;
+    float *t_input = calloc(t_intput_size, sizeof(float));
+    //char *
+    *t_bit_input = calloc(t_bit_input_size, sizeof(char));
+
+    //printf("\n bit_input_size = %d, n = %d, k = %d, ldb = %d \n", bit_input_size, n, k, n);
+    //printf("\n t_bit_input_size = %d, k = %d, n = %d, new_ldb = %d \n", t_bit_input_size, k, n, new_ldb);
+
+    //printf("\n align_weights_size = %d, k = %d, m = %d, lda = %d \n", align_weights_size, k, m, k);
+    //printf("\n align_bit_weights_size = %d, k = %d, m = %d, new_lda = %d \n", align_bit_weights_size, k, m, new_ldb);
+
+    // transpose and align B
+    int i, j;
+    //#pragma omp parallel for
+    /*
+    for (i = 0; i < n; ++i) {
+        for (j = 0; j < k; ++j) {
+            t_input[i*new_ldb + j] = b[j*n + i];
+        }
+    }*/
+    //transpose_block_SSE4x4(float *A, float *B, const int n, const int m, const int lda, const int ldb, const int block_size)
+
+    //transpose_block(b, t_input, k, n, n, new_ldb, 16);
+
+    int blocksize = 1;
+    int mod_k = 1, mod_n = 1;
+    for (i = 2; i < 256; i *= 2)
+        if (k % i == 0) mod_k = i;
+
+    for (i = 2; i < 256; i *= 2)
+        if (n % i == 0) mod_n = i;
+
+    blocksize = (mod_k < mod_n) ? mod_k : mod_n;
+
+    transpose_block_SSE4x4(b, t_input, k, n, n, new_ldb, blocksize);
+
+    //transpose_block(b, t_input, k, n, n, new_ldb, blocksize);
+    //printf("\n blocksize = %d \n", blocksize);
+
+    float_to_bit(t_input, *t_bit_input, t_intput_size);
+    free(t_input);
+
+    return t_intput_size;
+}
+
+
 void forward_convolutional_layer(convolutional_layer l, network_state state)
 {
     int out_h = convolutional_out_height(l);
@@ -652,8 +702,9 @@
     u++;
 
     for(i = 0; i < l.batch; ++i){
-        im2col_cpu(state.input, l.c, l.h, l.w,
-                l.size, l.stride, l.pad, b);
+        //im2col_cpu(state.input, l.c, l.h, l.w, l.size, l.stride, l.pad, b);
+        im2col_cpu_custom(state.input, l.c, l.h, l.w, l.size, l.stride, l.pad, b);
+
         //gemm(0,0,m,n,k,1,a,k,b,n,1,c,n);
         //gemm_nn_custom(m, n, k, 1, a, k, b, n, c, n);
         if (l.xnor) {
@@ -683,8 +734,8 @@
 
             // transpose B from NxK to KxN (x-axis (ldb = l.size*l.size*l.c) - should be multiple of 8 bits)
             {
+                /*
                 size_t ldb_align = 256;// 8;
-                if (k > 4096)ldb_align = 4096;
 
                 size_t new_ldb = k + (ldb_align - k%ldb_align); // (k / 8 + 1) * 8;
                 size_t t_intput_size = new_ldb * n;
@@ -709,6 +760,8 @@
                 }
                 float_to_bit(t_input, t_bit_input, t_intput_size);
 
+
+
                 if (!l.align_bit_weights)
                 {
                     size_t align_weights_size = new_ldb * m;
@@ -729,12 +782,17 @@
 
                     free(align_weights);
                 }
+                */
+                size_t ldb_align = 256; // 256 bit for AVX2
+                size_t new_ldb = k + (ldb_align - k%ldb_align);
+                char *t_bit_input = NULL;
+                size_t t_intput_size = binary_transpose_align_input(k, n, b, &t_bit_input, ldb_align);
 
                 gemm_nn_custom_bin_mean_transposed(m, n, k, 1, l.align_bit_weights, new_ldb, t_bit_input, new_ldb, c, n, l.mean_arr);
 
                 //gemm_nn_custom_bin_mean_transposed(m, n, k, 1, bit_weights, k, t_bit_input, new_ldb, c, n, mean_arr);
 
-                free(t_input);
+                //free(t_input);
                 free(t_bit_input);
 
                 //free(align_bit_weights);
@@ -771,7 +829,9 @@
     }
     add_bias(l.output, l.biases, l.batch, l.n, out_h*out_w);
 
-    activate_array(l.output, m*n*l.batch, l.activation);
+    //activate_array(l.output, m*n*l.batch, l.activation);
+    activate_array_cpu_custom(l.output, m*n*l.batch, l.activation);
+
     if(l.binary || l.xnor) swap_binary(&l);
 }
 
diff --git a/src/convolutional_layer.h b/src/convolutional_layer.h
index dd79c48..b804afb 100644
--- a/src/convolutional_layer.h
+++ b/src/convolutional_layer.h
@@ -35,7 +35,7 @@
 void swap_binary(convolutional_layer *l);
 void binarize_weights2(float *weights, int n, int size, char *binary, float *scales);
 
-void binary_transpose_align_weights(convolutional_layer *l, size_t ldb_align);
+void binary_align_weights(convolutional_layer *l, size_t ldb_align);
 
 void backward_convolutional_layer(convolutional_layer layer, network_state state);
 
diff --git a/src/gemm.c b/src/gemm.c
index 63163a5..478e966 100644
--- a/src/gemm.c
+++ b/src/gemm.c
@@ -1,5 +1,6 @@
 #include "gemm.h"
 #include "utils.h"
+#include "im2col.h"
 #include "cuda.h"
 #include <stdlib.h>
 #include <stdio.h>
@@ -426,7 +427,7 @@
 
 // http://graphics.stanford.edu/~seander/bithacks.html
 // https://stackoverflow.com/questions/17354971/fast-counting-the-number-of-set-bits-in-m128i-register
-
+// https://arxiv.org/pdf/1611.07612.pdf
 
 static inline int popcnt128(__m128i n) {
     const __m128i n_hi = _mm_unpackhi_epi64(n, n);
@@ -458,133 +459,61 @@
     return _mm256_sad_epu8(total, _mm256_setzero_si256());
 }

 static inline int popcnt256_custom(__m256i n) {
-    return _mm_popcnt_u64(n.m256i_i64[0]) +
-        _mm_popcnt_u64(n.m256i_i64[1]) +
-        _mm_popcnt_u64(n.m256i_i64[2]) +
-        _mm_popcnt_u64(n.m256i_i64[3]);
+    __m256i val = count256(n);
+
+    return val.m256i_i64[0] +
+    val.m256i_i64[1] +
+    val.m256i_i64[2] +
+    val.m256i_i64[3];
 }
 
-static inline void CSA(__m256i * h, __m256i * l, __m256i a, __m256i b, __m256i c)
-{
-    __m256i u = _mm256_xor_si256(a, b);
-    *h = _mm256_or_si256(_mm256_and_si256(a, b), _mm256_and_si256(u, c));
-    *l = _mm256_xor_si256(u, c);
-}
-
-static inline __m256i xnor256(__m256i a_bit256, __m256i b_bit256) {
-    __m256i all_1 = _mm256_set1_epi8(255);
-    __m256i xor256 = _mm256_xor_si256(a_bit256, b_bit256);
-    __m256i c_bit256 = _mm256_andnot_si256(xor256, all_1);
-
-    return c_bit256;
-
-}
-
-// 2 x faster than popcnt: https://arxiv.org/pdf/1611.07612.pdf
-// step = 16*256/8 = 512 bytes = 4096 bit (ldb, lda, bit_step, align - all should be aligned by 4096 bit)
-static inline uint64_t avx_hs_custom(__m256i * A, __m256i * B, uint64_t size) {
-    __m256i total = _mm256_setzero_si256();
-    __m256i ones = _mm256_setzero_si256();
-    __m256i twos = _mm256_setzero_si256();
-    __m256i fours = _mm256_setzero_si256();
-    __m256i eights = _mm256_setzero_si256();
-    __m256i sixteens = _mm256_setzero_si256();
-    __m256i twosA, twosB, foursA, foursB, eightsA, eightsB;
-
-    for (uint64_t i = 0; i < size; i += 16) {
-        //CSA(&twosA, &ones, ones, d[i], d[i + 1]);
-        CSA(&twosA, &ones, ones, xnor256(A[i], B[i]), xnor256(A[i + 1], B[i + 1]));
-        CSA(&twosB, &ones, ones, xnor256(A[i + 2], B[i + 2]), xnor256(A[i + 3], B[i + 3]));
-        CSA(&foursA, &twos, twos, twosA, twosB);
-        CSA(&twosA, &ones, ones, xnor256(A[i + 4], B[i + 4]), xnor256(A[i + 5], B[i + 5]));
-        CSA(&twosB, &ones, ones, xnor256(A[i + 6], B[i + 6]), xnor256(A[i + 7], B[i + 7]));
-        CSA(&foursB, &twos, twos, twosA, twosB);
-        CSA(&eightsA, &fours, fours, foursA, foursB);
-        CSA(&twosA, &ones, ones, xnor256(A[i + 8], B[i + 8]), xnor256(A[i + 9], B[i + 9]));
-        CSA(&twosB, &ones, ones, xnor256(A[i + 10], B[i + 10]), xnor256(A[i + 11], B[i + 11]));
-        CSA(&foursA, &twos, twos, twosA, twosB);
-        CSA(&twosA, &ones, ones, xnor256(A[i + 12], B[i + 12]), xnor256(A[i + 13], B[i + 13]));
-        CSA(&twosB, &ones, ones, xnor256(A[i + 14], B[i + 14]), xnor256(A[i + 15], B[i + 15]));
-        CSA(&foursB, &twos, twos, twosA, twosB);
-        CSA(&eightsB, &fours, fours, foursA, foursB);
-        CSA(&sixteens, &eights, eights, eightsA, eightsB);
-
-        total = _mm256_add_epi64(total, count256(sixteens));
-    }
-    total = _mm256_slli_epi64(total, 4);
-    total = _mm256_add_epi64(total,
-        _mm256_slli_epi64(count256(eights), 3));
-    total = _mm256_add_epi64(total,
-        _mm256_slli_epi64(count256(fours), 2));
-    total = _mm256_add_epi64(total,
-        _mm256_slli_epi64(count256(twos), 1));
-    total = _mm256_add_epi64(total, count256(ones));
-
-    return total.m256i_i64[0] +
-            total.m256i_i64[1] +
-            total.m256i_i64[2] +
-            total.m256i_i64[3];
-
-    //return _mm256_extract_epi64(total, 0)
-    //    + _mm256_extract_epi64(total, 1)
-    //    + _mm256_extract_epi64(total, 2)
-    //    + _mm256_extract_epi64(total, 3);
-}

-
 void gemm_nn_custom_bin_mean_transposed(int M, int N, int K, float ALPHA_UNUSED,
     unsigned char *A, int lda,
     unsigned char *B, int ldb,
     float *C, int ldc, float *mean_arr)
 {
-    __m256i all_1 = _mm256_set1_epi8(255);
-    int i, j, k;
+    int i;
 
-    //printf("\n M = %d, N = %d, K = %d, ldb = %d, M*ldb/8 = %d, N*ldb/8= %d \n", M, N, K, ldb, M*ldb/8, N*ldb/8);
-    //if (K > 4096)  printf("!!!avx_hs!!! \n\n");
+    static int max_num_threads = 0;
+    if (max_num_threads == 0) {
+        max_num_threads = omp_get_max_threads();
+        omp_set_num_threads(max_num_threads / 2);
+    }
 
     #pragma omp parallel for
-    for (i = 0; i < M; ++i) {   // l.n - filters [16 - 55 - 1024]
+    for (i = 0; i < M; ++i)
+    {   // l.n - filters [16 - 55 - 1024]
         float mean_val = mean_arr[i];
+        int j, k;
+        __m256i all_1 = _mm256_set1_epi8(255);
 
         for (j = 0; j < N; ++j) { // out_h*out_w - one channel output size [169 - 173056]
             int count = 0;
             const int bit_step = 256;
+            __m256i count_sum = _mm256_set1_epi8(0);
 
+            for (k = 0; k < K; k += bit_step) {   // l.size*l.size*l.c - one filter size [27 - 9216]
+                __m256i a_bit256 = _mm256_loadu_si256((__m256i *)(A + (i*lda + k) / 8));
+                __m256i b_bit256 = _mm256_loadu_si256((__m256i *)(B + (j*ldb + k) / 8));
+                __m256i xor256 = _mm256_xor_si256(a_bit256, b_bit256);  // xnor = not(xor(a,b))
+                __m256i c_bit256 = _mm256_andnot_si256(xor256, all_1);  // can be optimized - we can do other NOT for wegihts once and do not do this NOT
 
-            int hs_count = 0;
-            if (K > 4096) {
-                hs_count = avx_hs_custom(A + (i*lda) / 8, B + (j*ldb) / 8, K / 256);
+                count_sum = _mm256_add_epi64(count256(c_bit256), count_sum);    //  Mula�s algorithm
 
-                int local_bit_step = 4096;
+                //count += popcnt256(c_bit256);
 
-                int f1 = (K % local_bit_step == 0) ? 0 : (local_bit_step - (K % local_bit_step));
-                hs_count = hs_count - f1;    // remove extra bits
-                count = hs_count;
+                //binary_int64_printf(c_bit64);
+                //printf(", count = %d \n\n", tmp_count);
             }
-            else {
-                for (k = 0; k < K; k += bit_step) {   // l.size*l.size*l.c - one filter size [27 - 9216]
 
-                    //__m128i a_bit128 = _mm_loadu_si128((__m128i *)(A + (i*lda + k) / 8));
-                    //__m128i b_bit128 = _mm_loadu_si128((__m128i *)(B + (j*ldb + k) / 8));
-                    //__m128i xor128 = _mm_xor_si128(a_bit128, b_bit128);
-                    //__m128i c_bit128 = _mm_andnot_si128(xor128, all_1);
-                    //int tmp_count = popcnt128(c_bit128);
+            // count of 1 bits
+            count = count_sum.m256i_i64[0] +
+                count_sum.m256i_i64[1] +
+                count_sum.m256i_i64[2] +
+                count_sum.m256i_i64[3];
 
-                    __m256i a_bit256 = _mm256_loadu_si256((__m256i *)(A + (i*lda + k) / 8));
-                    __m256i b_bit256 = _mm256_loadu_si256((__m256i *)(B + (j*ldb + k) / 8));
-                    __m256i xor256 = _mm256_xor_si256(a_bit256, b_bit256);
-                    __m256i c_bit256 = _mm256_andnot_si256(xor256, all_1); //we can do NOT for wegihts once and do not do this NOT
-                    int tmp_count = popcnt256(c_bit256);
-                    //int tmp_count = popcnt256_custom(c_bit256);
-                    count += tmp_count;
-
-                    //binary_int64_printf(c_bit64);
-                    //printf(", count = %d \n\n", tmp_count);
-                }
-
-                int f1 = (K % bit_step == 0) ? 0 : (bit_step - (K % bit_step));
-                count = count - f1;    // remove extra bits
-           }
+            int f1 = (K % bit_step == 0) ? 0 : (bit_step - (K % bit_step));
+            count = count - f1;    // remove extra bits (from empty space for align only)
 
             C[i*ldc + j] = (2 * count - K) * mean_val;
         }
@@ -592,6 +521,142 @@
 }
 
 
+static inline float im2col_get_pixel(float *im, int height, int width, int channels,
+    int row, int col, int channel, int pad)
+{
+    row -= pad;
+    col -= pad;
+
+    if (row < 0 || col < 0 ||
+        row >= height || col >= width) return 0;
+    return im[col + width*(row + height*channel)];
+}
+
+//From Berkeley Vision's Caffe!
+//https://github.com/BVLC/caffe/blob/master/LICENSE
+void im2col_cpu_custom(float* data_im,
+    int channels, int height, int width,
+    int ksize, int stride, int pad, float* data_col)
+{
+
+    int c, h, w;
+    int height_col = (height + 2 * pad - ksize) / stride + 1;
+    int width_col = (width + 2 * pad - ksize) / stride + 1;
+    int channels_col = channels * ksize * ksize;
+
+    // optimized version
+    if (height_col == height && width_col == width && stride == 1 && pad == 1)
+    {
+        #pragma omp parallel for
+        for (c = 0; c < channels_col; ++c) {
+            int w_offset = c % ksize;
+            int h_offset = (c / ksize) % ksize;
+            int c_im = c / ksize / ksize;
+            for (h = pad; h < height_col-pad; ++h) {
+                for (w = pad; w < width_col-pad-8; w += 8) {
+                    int im_row = h_offset + h - pad;
+                    int im_col = w_offset + w - pad;
+                    int col_index = (c * height_col + h) * width_col + w;
+
+                    //data_col[col_index] = data_im[im_col + width*(im_row + height*c_im)];
+                    __m256 src256 = _mm256_loadu_ps((__m256i *)(&data_im[im_col + width*(im_row + height*c_im)]));
+                    _mm256_storeu_ps(&data_col[col_index], src256);
+                }
+
+                for (; w < width_col - pad; ++w) {
+                    int im_row = h_offset + h - pad;
+                    int im_col = w_offset + w - pad;
+                    int col_index = (c * height_col + h) * width_col + w;
+
+                    data_col[col_index] = data_im[im_col + width*(im_row + height*c_im)];
+                }
+            }
+
+            {
+                w = 0;
+                for (h = 0; h < height_col; ++h) {
+                    int im_row = h_offset + h;
+                    int im_col = w_offset + w;
+                    int col_index = (c * height_col + h) * width_col + w;
+                    data_col[col_index] = im2col_get_pixel(data_im, height, width, channels,
+                        im_row, im_col, c_im, pad);
+                }
+            }
+
+            {
+                w = width_col-1;
+                for (h = 0; h < height_col; ++h) {
+                    int im_row = h_offset + h;
+                    int im_col = w_offset + w;
+                    int col_index = (c * height_col + h) * width_col + w;
+                    data_col[col_index] = im2col_get_pixel(data_im, height, width, channels,
+                        im_row, im_col, c_im, pad);
+                }
+            }
+
+            {
+                h = 0;
+                for (w = 0; w < width_col; ++w) {
+                    int im_row = h_offset + h;
+                    int im_col = w_offset + w;
+                    int col_index = (c * height_col + h) * width_col + w;
+                    data_col[col_index] = im2col_get_pixel(data_im, height, width, channels,
+                            im_row, im_col, c_im, pad);
+                }
+            }
+
+            {
+                h = height_col-1;
+                for (w = 0; w < width_col; ++w) {
+                    int im_row = h_offset + h;
+                    int im_col = w_offset + w;
+                    int col_index = (c * height_col + h) * width_col + w;
+                    data_col[col_index] = im2col_get_pixel(data_im, height, width, channels,
+                        im_row, im_col, c_im, pad);
+                }
+            }
+        }
+
+    }
+    else {
+        //printf("\n Error: is no non-optimized version \n");
+        im2col_cpu(data_im, channels, height, width, ksize, stride, pad, data_col);
+    }
+}
+
+void activate_array_cpu_custom(float *x, const int n, const ACTIVATION a)
+{
+    int i;
+    if (a == LINEAR)
+    {}
+    else if (a == LEAKY)
+    {
+        __m256i all256_sing1 = _mm256_set_epi32(0x80000000, 0x80000000, 0x80000000, 0x80000000, 0x80000000, 0x80000000, 0x80000000, 0x80000000);
+        __m256 all256_01 = _mm256_set1_ps(0.1F);
+
+        for (i = 0; i < n; i += 8) {
+            //x[i] = (x[i]>0) ? x[i] : .1*x[i];
+
+            __m256 src256 = _mm256_loadu_ps((__m256 *)(&x[i]));
+            __m256 mult256 = _mm256_mul_ps((src256), all256_01); // mult * 0.1
+
+            __m256i sign256 = _mm256_and_si256(_mm256_castps_si256(src256), all256_sing1); // check sign in 8 x 32-bit floats
+
+            __m256 result256 = _mm256_blendv_ps(src256, mult256, _mm256_castsi256_ps(sign256)); // (sign>0) ? src : mult;
+            _mm256_storeu_ps((__m256 *)(&x[i]), result256);
+        }
+
+        for (; i < n; ++i) {
+            x[i] = (x[i]>0) ? x[i] : .1*x[i];
+        }
+    }
+    else {
+        for (i = 0; i < n; ++i) {
+            x[i] = activate(x[i], a);
+        }
+    }
+}
+
 void float_to_bit(float *src, unsigned char *dst, size_t size)
 {
     size_t dst_size = size / 8 + 1;
@@ -612,6 +677,56 @@
     }
 }
 
+static inline void transpose4x4_SSE(float *A, float *B, const int lda, const int ldb)
+{
+    __m128 row1 = _mm_load_ps(&A[0 * lda]);
+    __m128 row2 = _mm_load_ps(&A[1 * lda]);
+    __m128 row3 = _mm_load_ps(&A[2 * lda]);
+    __m128 row4 = _mm_load_ps(&A[3 * lda]);
+    _MM_TRANSPOSE4_PS(row1, row2, row3, row4);
+    _mm_store_ps(&B[0 * ldb], row1);
+    _mm_store_ps(&B[1 * ldb], row2);
+    _mm_store_ps(&B[2 * ldb], row3);
+    _mm_store_ps(&B[3 * ldb], row4);
+}
+
+void transpose_block_SSE4x4(float *A, float *B, const int n, const int m,
+    const int lda, const int ldb, const int block_size)
+{
+    int i;
+    if (block_size % 4 == 0) {
+        #pragma omp parallel for
+        for (i = 0; i < n; i += block_size) {
+            int j, i2, j2;
+            for (j = 0; j < m; j += block_size) {
+                int max_i2 = i + block_size < n ? i + block_size : n;
+                int max_j2 = j + block_size < m ? j + block_size : m;
+                for (i2 = i; i2 < max_i2; i2 += 4) {
+                    for (j2 = j; j2 < max_j2; j2 += 4) {
+                        transpose4x4_SSE(&A[i2*lda + j2], &B[j2*ldb + i2], lda, ldb);
+                    }
+                }
+            }
+        }
+    }
+    else {
+        #pragma omp parallel for
+        for (i = 0; i < n; i += block_size) {
+            int j, i2, j2;
+            for (j = 0; j < m; j += block_size) {
+                int max_i2 = i + block_size < n ? i + block_size : n;
+                int max_j2 = j + block_size < m ? j + block_size : m;
+                for (i2 = i; i2 < max_i2; ++i2) {
+                    for (j2 = j; j2 < max_j2; ++j2) {
+                        B[j2*ldb + i2] = A[i2*lda + j2];
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 #else
 
 void gemm_nn(int M, int N, int K, float ALPHA,
@@ -666,6 +781,115 @@
     }
 }
 
+//From Berkeley Vision's Caffe!
+//https://github.com/BVLC/caffe/blob/master/LICENSE
+void im2col_cpu_custom(float* data_im,
+    int channels, int height, int width,
+    int ksize, int stride, int pad, float* data_col)
+{
+
+    int c, h, w;
+    int height_col = (height + 2 * pad - ksize) / stride + 1;
+    int width_col = (width + 2 * pad - ksize) / stride + 1;
+    int channels_col = channels * ksize * ksize;
+
+    // optimized version
+    if (height_col == height && width_col == width && stride == 1 && pad == 1)
+    {
+        #pragma omp parallel for
+        for (c = 0; c < channels_col; ++c) {
+            int w_offset = c % ksize;
+            int h_offset = (c / ksize) % ksize;
+            int c_im = c / ksize / ksize;
+            for (h = pad; h < height_col - pad; ++h) {
+                for (w = pad; w < width_col - pad; ++w) {
+                    int im_row = h_offset + h - pad;
+                    int im_col = w_offset + w - pad;
+                    int col_index = (c * height_col + h) * width_col + w;
+
+                    data_col[col_index] = data_im[im_col + width*(im_row + height*c_im)];
+    }
+
+                for (; w < width_col - pad; ++w) {
+                    int im_row = h_offset + h - pad;
+                    int im_col = w_offset + w - pad;
+                    int col_index = (c * height_col + h) * width_col + w;
+
+                    data_col[col_index] = data_im[im_col + width*(im_row + height*c_im)];
+                }
+}
+
+            {
+                w = 0;
+                for (h = 0; h < height_col; ++h) {
+                    int im_row = h_offset + h;
+                    int im_col = w_offset + w;
+                    int col_index = (c * height_col + h) * width_col + w;
+                    data_col[col_index] = im2col_get_pixel(data_im, height, width, channels,
+                        im_row, im_col, c_im, pad);
+                }
+            }
+
+            {
+                w = width_col - 1;
+                for (h = 0; h < height_col; ++h) {
+                    int im_row = h_offset + h;
+                    int im_col = w_offset + w;
+                    int col_index = (c * height_col + h) * width_col + w;
+                    data_col[col_index] = im2col_get_pixel(data_im, height, width, channels,
+                        im_row, im_col, c_im, pad);
+                }
+            }
+
+            {
+                h = 0;
+                for (w = 0; w < width_col; ++w) {
+                    int im_row = h_offset + h;
+                    int im_col = w_offset + w;
+                    int col_index = (c * height_col + h) * width_col + w;
+                    data_col[col_index] = im2col_get_pixel(data_im, height, width, channels,
+                        im_row, im_col, c_im, pad);
+                }
+            }
+
+            {
+                h = height_col - 1;
+                for (w = 0; w < width_col; ++w) {
+                    int im_row = h_offset + h;
+                    int im_col = w_offset + w;
+                    int col_index = (c * height_col + h) * width_col + w;
+                    data_col[col_index] = im2col_get_pixel(data_im, height, width, channels,
+                        im_row, im_col, c_im, pad);
+                }
+            }
+        }
+
+    }
+    else {
+        //printf("\n Error: is no non-optimized version \n");
+        im2col_cpu(data_im, channels, height, width, ksize, stride, pad, data_col);
+    }
+}
+
+void activate_array_cpu_custom(float *x, const int n, const ACTIVATION a)
+{
+    int i;
+    if (a == LINEAR)
+    {
+    }
+    else if (a == LEAKY)
+    {
+        for (i = 0; i < n; ++i) {
+            x[i] = (x[i]>0) ? x[i] : .1*x[i];
+        }
+    }
+    else {
+        for (i = 0; i < n; ++i) {
+            x[i] = activate(x[i], a);
+        }
+    }
+}
+
 void float_to_bit(float *src, unsigned char *dst, size_t size)
 {
     size_t dst_size = size / 8 + 1;
@@ -695,6 +919,36 @@
     }
     free(byte_arr);
 }
+
+static inline void transpose_scalar_block(float *A, float *B, const int lda, const int ldb, const int block_size)
+{
+    int i, j;
+    //#pragma omp parallel for
+    for (i = 0; i<block_size; i++) {
+        for (j = 0; j<block_size; j++) {
+            B[j*ldb + i] = A[i*lda + j];
+        }
+    }
+}
+
+void transpose_block_SSE4x4(float *A, float *B, const int n, const int m,
+    const int lda, const int ldb, const int block_size)
+{
+    int i;
+    #pragma omp parallel for
+    for (i = 0; i < n; i += block_size) {
+        int j, i2, j2;
+        for (j = 0; j < m; j += block_size) {
+            int max_i2 = i + block_size < n ? i + block_size : n;
+            int max_j2 = j + block_size < m ? j + block_size : m;
+            for (i2 = i; i2 < max_i2; ++i2) {
+                for (j2 = j; j2 < max_j2; ++j2) {
+                    B[j2*ldb + i2] = A[i2*lda + j2];
+                }
+                }
+            }
+        }
+}
 #endif    // __x86_64
 
 void gemm_nt(int M, int N, int K, float ALPHA,
diff --git a/src/gemm.h b/src/gemm.h
index 8664d8e..97fa09c 100644
--- a/src/gemm.h
+++ b/src/gemm.h
@@ -1,5 +1,6 @@
 #ifndef GEMM_H
 #define GEMM_H
+#include "activations.h"
 
 static inline void set_bit(unsigned char *const dst, size_t index) {
     size_t dst_i = index / 8;
@@ -16,17 +17,19 @@
 
 void float_to_bit(float *src, unsigned char *dst, size_t size);
 
+void transpose_block_SSE4x4(float *A, float *B, const int n, const int m,
+    const int lda, const int ldb, const int block_size);
+
 void gemm_nn_custom_bin_mean_transposed(int M, int N, int K, float ALPHA_UNUSED,
     unsigned char *A, int lda,
     unsigned char *B, int ldb,
     float *C, int ldc, float *mean_arr);
 
+void im2col_cpu_custom(float* data_im,
+    int channels, int height, int width,
+    int ksize, int stride, int pad, float* data_col);
 
-//void gemm_nn_custom_bin_mean(int M, int N, int K, float ALPHA_UNUSED,
-    //unsigned char *A, int lda,
-    //unsigned char *B, int ldb,
-    //float *C, int ldc, float *mean_arr)
-
+void activate_array_cpu_custom(float *x, const int n, const ACTIVATION a);
 
 
 void gemm_bin(int M, int N, int K, float ALPHA,
diff --git a/src/network.c b/src/network.c
index a62f6d0..345ce68 100644
--- a/src/network.c
+++ b/src/network.c
@@ -862,9 +862,8 @@
             if (l->xnor) {
                 //printf("\n %d \n", j);
                 size_t ldb_align = 256; // 256bit for AVX2
-                if (l->size*l->size*l->c > 4096) ldb_align = 4096;
 
-                binary_transpose_align_weights(l, ldb_align);
+                binary_align_weights(l, ldb_align);
             }
         }
     }

--
Gitblit v1.10.0