Squashed 'third_party/blasfeo/' content from commit 2a828ca
Change-Id: If1c3caa4799b2d4eb287ef83fa17043587ef07a3
git-subtree-dir: third_party/blasfeo
git-subtree-split: 2a828ca5442108c4c58e4b42b061a0469043f6ea
diff --git a/kernel/avx/kernel_sgemm_diag_lib8.c b/kernel/avx/kernel_sgemm_diag_lib8.c
new file mode 100644
index 0000000..63183b2
--- /dev/null
+++ b/kernel/avx/kernel_sgemm_diag_lib8.c
@@ -0,0 +1,480 @@
+/**************************************************************************************************
+* *
+* This file is part of BLASFEO. *
+* *
+* BLASFEO -- BLAS For Embedded Optimization. *
+* Copyright (C) 2016-2017 by Gianluca Frison. *
+* Developed at IMTEK (University of Freiburg) under the supervision of Moritz Diehl. *
+* All rights reserved. *
+* *
+* HPMPC is free software; you can redistribute it and/or *
+* modify it under the terms of the GNU Lesser General Public *
+* License as published by the Free Software Foundation; either *
+* version 2.1 of the License, or (at your option) any later version. *
+* *
+* HPMPC is distributed in the hope that it will be useful, *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
+* See the GNU Lesser General Public License for more details. *
+* *
+* You should have received a copy of the GNU Lesser General Public *
+* License along with HPMPC; if not, write to the Free Software *
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
+* *
+* Author: Gianluca Frison, giaf (at) dtu.dk *
+* gianluca.frison (at) imtek.uni-freiburg.de *
+* *
+**************************************************************************************************/
+
+#include <mmintrin.h>
+#include <xmmintrin.h> // SSE
+#include <emmintrin.h> // SSE2
+#include <pmmintrin.h> // SSE3
+#include <smmintrin.h> // SSE4
+#include <immintrin.h> // AVX
+
+
+
+// B is the diagonal of a matrix, beta==0.0 case
+void kernel_sgemm_diag_right_4_a0_lib4(int kmax, float *alpha, float *A, int sda, float *B, float *D, int sdd)
+ {
+
+ if(kmax<=0)
+ return;
+
+ const int bs = 8;
+
+ int k;
+
+ __m256
+ alpha0,
+ mask_f,
+ sign,
+ a_00,
+ b_00, b_11, b_22, b_33,
+ d_00, d_01, d_02, d_03;
+
+ __m256i
+ mask_i;
+
+ alpha0 = _mm256_broadcast_ss( alpha );
+
+ b_00 = _mm256_broadcast_ss( &B[0] );
+ b_00 = _mm256_mul_ps( b_00, alpha0 );
+ b_11 = _mm256_broadcast_ss( &B[1] );
+ b_11 = _mm256_mul_ps( b_11, alpha0 );
+ b_22 = _mm256_broadcast_ss( &B[2] );
+ b_22 = _mm256_mul_ps( b_22, alpha0 );
+ b_33 = _mm256_broadcast_ss( &B[3] );
+ b_33 = _mm256_mul_ps( b_33, alpha0 );
+
+ for(k=0; k<kmax-7; k+=8)
+ {
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+ a_00 = _mm256_load_ps( &A[8] );
+ d_01 = _mm256_mul_ps( a_00, b_11 );
+ a_00 = _mm256_load_ps( &A[16] );
+ d_02 = _mm256_mul_ps( a_00, b_22 );
+ a_00 = _mm256_load_ps( &A[24] );
+ d_03 = _mm256_mul_ps( a_00, b_33 );
+
+ _mm256_store_ps( &D[0], d_00 );
+ _mm256_store_ps( &D[8], d_01 );
+ _mm256_store_ps( &D[16], d_02 );
+ _mm256_store_ps( &D[24], d_03 );
+
+ A += 8*sda;
+ D += 8*sdd;
+
+ }
+ if(k<kmax)
+ {
+
+ const float mask_f[] = {0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5};
+ float m_f = kmax-k;
+
+ mask_i = _mm256_castps_si256( _mm256_sub_ps( _mm256_loadu_ps( mask_f ), _mm256_broadcast_ss( &m_f ) ) );
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+ a_00 = _mm256_load_ps( &A[8] );
+ d_01 = _mm256_mul_ps( a_00, b_11 );
+ a_00 = _mm256_load_ps( &A[16] );
+ d_02 = _mm256_mul_ps( a_00, b_22 );
+ a_00 = _mm256_load_ps( &A[24] );
+ d_03 = _mm256_mul_ps( a_00, b_33 );
+
+ _mm256_maskstore_ps( &D[0], mask_i, d_00 );
+ _mm256_maskstore_ps( &D[8], mask_i, d_01 );
+ _mm256_maskstore_ps( &D[16], mask_i, d_02 );
+ _mm256_maskstore_ps( &D[24], mask_i, d_03 );
+
+ }
+
+ }
+
+
+
+// B is the diagonal of a matrix
+void kernel_sgemm_diag_right_4_lib4(int kmax, float *alpha, float *A, int sda, float *B, float *beta, float *C, int sdc, float *D, int sdd)
+ {
+
+ if(kmax<=0)
+ return;
+
+ const int bs = 8;
+
+ int k;
+
+ __m256
+ alpha0, beta0,
+ mask_f,
+ sign,
+ a_00,
+ b_00, b_11, b_22, b_33,
+ c_00,
+ d_00, d_01, d_02, d_03;
+
+ __m256i
+ mask_i;
+
+ alpha0 = _mm256_broadcast_ss( alpha );
+ beta0 = _mm256_broadcast_ss( beta );
+
+ b_00 = _mm256_broadcast_ss( &B[0] );
+ b_00 = _mm256_mul_ps( b_00, alpha0 );
+ b_11 = _mm256_broadcast_ss( &B[1] );
+ b_11 = _mm256_mul_ps( b_11, alpha0 );
+ b_22 = _mm256_broadcast_ss( &B[2] );
+ b_22 = _mm256_mul_ps( b_22, alpha0 );
+ b_33 = _mm256_broadcast_ss( &B[3] );
+ b_33 = _mm256_mul_ps( b_33, alpha0 );
+
+ for(k=0; k<kmax-7; k+=8)
+ {
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+ a_00 = _mm256_load_ps( &A[8] );
+ d_01 = _mm256_mul_ps( a_00, b_11 );
+ a_00 = _mm256_load_ps( &A[16] );
+ d_02 = _mm256_mul_ps( a_00, b_22 );
+ a_00 = _mm256_load_ps( &A[24] );
+ d_03 = _mm256_mul_ps( a_00, b_33 );
+
+ c_00 = _mm256_load_ps( &C[0] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_00 = _mm256_add_ps( c_00, d_00 );
+ c_00 = _mm256_load_ps( &C[8] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_01 = _mm256_add_ps( c_00, d_01 );
+ c_00 = _mm256_load_ps( &C[16] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_02 = _mm256_add_ps( c_00, d_02 );
+ c_00 = _mm256_load_ps( &C[24] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_03 = _mm256_add_ps( c_00, d_03 );
+
+ _mm256_store_ps( &D[0], d_00 );
+ _mm256_store_ps( &D[8], d_01 );
+ _mm256_store_ps( &D[16], d_02 );
+ _mm256_store_ps( &D[24], d_03 );
+
+ A += 8*sda;
+ C += 8*sdc;
+ D += 8*sdd;
+
+ }
+ if(k<kmax)
+ {
+
+ const float mask_f[] = {0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5};
+ float m_f = kmax-k;
+
+ mask_i = _mm256_castps_si256( _mm256_sub_ps( _mm256_loadu_ps( mask_f ), _mm256_broadcast_ss( &m_f ) ) );
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+ a_00 = _mm256_load_ps( &A[8] );
+ d_01 = _mm256_mul_ps( a_00, b_11 );
+ a_00 = _mm256_load_ps( &A[16] );
+ d_02 = _mm256_mul_ps( a_00, b_22 );
+ a_00 = _mm256_load_ps( &A[24] );
+ d_03 = _mm256_mul_ps( a_00, b_33 );
+
+ c_00 = _mm256_load_ps( &C[0] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_00 = _mm256_add_ps( c_00, d_00 );
+ c_00 = _mm256_load_ps( &C[8] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_01 = _mm256_add_ps( c_00, d_01 );
+ c_00 = _mm256_load_ps( &C[16] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_02 = _mm256_add_ps( c_00, d_02 );
+ c_00 = _mm256_load_ps( &C[24] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_03 = _mm256_add_ps( c_00, d_03 );
+
+ _mm256_maskstore_ps( &D[0], mask_i, d_00 );
+ _mm256_maskstore_ps( &D[8], mask_i, d_01 );
+ _mm256_maskstore_ps( &D[16], mask_i, d_02 );
+ _mm256_maskstore_ps( &D[24], mask_i, d_03 );
+
+ }
+
+ }
+
+
+
+// B is the diagonal of a matrix
+void kernel_sgemm_diag_right_3_lib4(int kmax, float *alpha, float *A, int sda, float *B, float *beta, float *C, int sdc, float *D, int sdd)
+ {
+
+ if(kmax<=0)
+ return;
+
+ const int bs = 8;
+
+ int k;
+
+ __m256
+ alpha0, beta0,
+ mask_f,
+ sign,
+ a_00,
+ b_00, b_11, b_22,
+ c_00,
+ d_00, d_01, d_02;
+
+ __m256i
+ mask_i;
+
+ alpha0 = _mm256_broadcast_ss( alpha );
+ beta0 = _mm256_broadcast_ss( beta );
+
+ b_00 = _mm256_broadcast_ss( &B[0] );
+ b_00 = _mm256_mul_ps( b_00, alpha0 );
+ b_11 = _mm256_broadcast_ss( &B[1] );
+ b_11 = _mm256_mul_ps( b_11, alpha0 );
+ b_22 = _mm256_broadcast_ss( &B[2] );
+ b_22 = _mm256_mul_ps( b_22, alpha0 );
+
+ for(k=0; k<kmax-7; k+=8)
+ {
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+ a_00 = _mm256_load_ps( &A[8] );
+ d_01 = _mm256_mul_ps( a_00, b_11 );
+ a_00 = _mm256_load_ps( &A[16] );
+ d_02 = _mm256_mul_ps( a_00, b_22 );
+
+ c_00 = _mm256_load_ps( &C[0] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_00 = _mm256_add_ps( c_00, d_00 );
+ c_00 = _mm256_load_ps( &C[8] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_01 = _mm256_add_ps( c_00, d_01 );
+ c_00 = _mm256_load_ps( &C[16] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_02 = _mm256_add_ps( c_00, d_02 );
+
+ _mm256_store_ps( &D[0], d_00 );
+ _mm256_store_ps( &D[8], d_01 );
+ _mm256_store_ps( &D[16], d_02 );
+
+ A += 8*sda;
+ C += 8*sdc;
+ D += 8*sdd;
+
+ }
+ if(k<kmax)
+ {
+
+ const float mask_f[] = {0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5};
+ float m_f = kmax-k;
+
+ mask_i = _mm256_castps_si256( _mm256_sub_ps( _mm256_loadu_ps( mask_f ), _mm256_broadcast_ss( &m_f ) ) );
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+ a_00 = _mm256_load_ps( &A[8] );
+ d_01 = _mm256_mul_ps( a_00, b_11 );
+ a_00 = _mm256_load_ps( &A[16] );
+ d_02 = _mm256_mul_ps( a_00, b_22 );
+
+ c_00 = _mm256_load_ps( &C[0] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_00 = _mm256_add_ps( c_00, d_00 );
+ c_00 = _mm256_load_ps( &C[8] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_01 = _mm256_add_ps( c_00, d_01 );
+ c_00 = _mm256_load_ps( &C[16] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_02 = _mm256_add_ps( c_00, d_02 );
+
+ _mm256_maskstore_ps( &D[0], mask_i, d_00 );
+ _mm256_maskstore_ps( &D[8], mask_i, d_01 );
+ _mm256_maskstore_ps( &D[16], mask_i, d_02 );
+
+ }
+
+ }
+
+
+
+// B is the diagonal of a matrix
+void kernel_sgemm_diag_right_2_lib4(int kmax, float *alpha, float *A, int sda, float *B, float *beta, float *C, int sdc, float *D, int sdd)
+ {
+
+ if(kmax<=0)
+ return;
+
+ const int bs = 4;
+
+ int k;
+
+ __m256
+ alpha0, beta0,
+ mask_f,
+ sign,
+ a_00,
+ b_00, b_11,
+ c_00,
+ d_00, d_01;
+
+ __m256i
+ mask_i;
+
+ alpha0 = _mm256_broadcast_ss( alpha );
+ beta0 = _mm256_broadcast_ss( beta );
+
+ b_00 = _mm256_broadcast_ss( &B[0] );
+ b_00 = _mm256_mul_ps( b_00, alpha0 );
+ b_11 = _mm256_broadcast_ss( &B[1] );
+ b_11 = _mm256_mul_ps( b_11, alpha0 );
+
+ for(k=0; k<kmax-7; k+=8)
+ {
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+ a_00 = _mm256_load_ps( &A[8] );
+ d_01 = _mm256_mul_ps( a_00, b_11 );
+
+ c_00 = _mm256_load_ps( &C[0] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_00 = _mm256_add_ps( c_00, d_00 );
+ c_00 = _mm256_load_ps( &C[8] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_01 = _mm256_add_ps( c_00, d_01 );
+
+ _mm256_store_ps( &D[0], d_00 );
+ _mm256_store_ps( &D[8], d_01 );
+
+ A += 8*sda;
+ C += 8*sdc;
+ D += 8*sdd;
+
+ }
+ if(k<kmax)
+ {
+
+ const float mask_f[] = {0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5};
+ float m_f = kmax-k;
+
+ mask_i = _mm256_castps_si256( _mm256_sub_ps( _mm256_loadu_ps( mask_f ), _mm256_broadcast_ss( &m_f ) ) );
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+ a_00 = _mm256_load_ps( &A[8] );
+ d_01 = _mm256_mul_ps( a_00, b_11 );
+
+ c_00 = _mm256_load_ps( &C[0] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_00 = _mm256_add_ps( c_00, d_00 );
+ c_00 = _mm256_load_ps( &C[8] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_01 = _mm256_add_ps( c_00, d_01 );
+
+ _mm256_maskstore_ps( &D[0], mask_i, d_00 );
+ _mm256_maskstore_ps( &D[8], mask_i, d_01 );
+
+ }
+
+ }
+
+
+
+// B is the diagonal of a matrix
+void kernel_sgemm_diag_right_1_lib4(int kmax, float *alpha, float *A, int sda, float *B, float *beta, float *C, int sdc, float *D, int sdd)
+ {
+
+ if(kmax<=0)
+ return;
+
+ const int bs = 4;
+
+ int k;
+
+ __m256
+ alpha0, beta0,
+ mask_f,
+ sign,
+ a_00,
+ b_00,
+ c_00,
+ d_00;
+
+ __m256i
+ mask_i;
+
+ alpha0 = _mm256_broadcast_ss( alpha );
+ beta0 = _mm256_broadcast_ss( beta );
+
+ b_00 = _mm256_broadcast_ss( &B[0] );
+ b_00 = _mm256_mul_ps( b_00, alpha0 );
+
+ for(k=0; k<kmax-7; k+=8)
+ {
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+
+ c_00 = _mm256_load_ps( &C[0] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_00 = _mm256_add_ps( c_00, d_00 );
+
+ _mm256_store_ps( &D[0], d_00 );
+
+ A += 8*sda;
+ C += 8*sdc;
+ D += 8*sdd;
+
+ }
+ if(k<kmax)
+ {
+
+ const float mask_f[] = {0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5};
+ float m_f = kmax-k;
+
+ mask_i = _mm256_castps_si256( _mm256_sub_ps( _mm256_loadu_ps( mask_f ), _mm256_broadcast_ss( &m_f ) ) );
+
+ a_00 = _mm256_load_ps( &A[0] );
+ d_00 = _mm256_mul_ps( a_00, b_00 );
+
+ c_00 = _mm256_load_ps( &C[0] );
+ c_00 = _mm256_mul_ps( c_00, beta0 );
+ d_00 = _mm256_add_ps( c_00, d_00 );
+
+ _mm256_maskstore_ps( &D[0], mask_i, d_00 );
+
+ }
+
+ }
+
+
+
+