#include "sn_intrinsic.h" #define _FMA_SELECT(i) _MM_SHUFFLE(i, i, i, i) // ans = left * right + ans // right is not transposed DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, sse, fma) { const __m128 a0 = PANEL_COL_VEC((*left), 0, 0); const __m128 a1 = PANEL_COL_VEC((*left), 1, 0); const __m128 a2 = PANEL_COL_VEC((*left), 2, 0); const __m128 a3 = PANEL_COL_VEC((*left), 3, 0); __m128 c0 = PANEL_COL_VEC((*ans), 0, 0); __m128 c1 = PANEL_COL_VEC((*ans), 1, 0); __m128 c2 = PANEL_COL_VEC((*ans), 2, 0); __m128 c3 = PANEL_COL_VEC((*ans), 3, 0); __m128 bcol; /* Column 0 */ bcol = PANEL_COL_VEC((*right), 0, 0); c0 = _mm_add_ps(c0, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0)))); c0 = _mm_add_ps(c0, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1)))); c0 = _mm_add_ps(c0, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2)))); c0 = _mm_add_ps(c0, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3)))); /* Column 1 */ bcol = PANEL_COL_VEC((*right), 1, 0); c1 = _mm_add_ps(c1, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0)))); c1 = _mm_add_ps(c1, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1)))); c1 = _mm_add_ps(c1, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2)))); c1 = _mm_add_ps(c1, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3)))); /* Column 2 */ bcol = PANEL_COL_VEC((*right), 2, 0); c2 = _mm_add_ps(c2, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0)))); c2 = _mm_add_ps(c2, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1)))); c2 = _mm_add_ps(c2, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2)))); c2 = _mm_add_ps(c2, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3)))); /* Column 3 */ bcol = PANEL_COL_VEC((*right), 3, 0); c3 = _mm_add_ps(c3, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0)))); c3 = _mm_add_ps(c3, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1)))); c3 = _mm_add_ps(c3, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2)))); c3 = _mm_add_ps(c3, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3)))); PANEL_COL_VEC((*ans), 0, 0) = c0; PANEL_COL_VEC((*ans), 1, 0) = c1; PANEL_COL_VEC((*ans), 2, 0) = c2; PANEL_COL_VEC((*ans), 3, 0) = c3; } // ans = left * right + ans // right is not transposed DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, avx_fma, fma) { static const __m256i idx0 = { 1, 1, 1, 1, 1, 1, 1, 1 }; for (int col = 0; col < 8; ++col) { register __m256 bcol = _mm256_load_ps((float *)&right->vec[col][0]); register __m256 ccol = _mm256_load_ps((float *)&ans->vec[col][0]); register __m256i idx = _mm256_setzero_si256(); register __m256i increament = _mm256_load_si256(&idx0); register __m256 a0 = _mm256_load_ps((float *)&left->vec[0][0]); register __m256 a1 = _mm256_load_ps((float *)&left->vec[1][0]); register __m256 a2 = _mm256_load_ps((float *)&left->vec[2][0]); register __m256 a3 = _mm256_load_ps((float *)&left->vec[3][0]); ccol = _mm256_fmadd_ps(a0, _mm256_permutevar8x32_ps(bcol, idx), ccol); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_fmadd_ps(a1, _mm256_permutevar8x32_ps(bcol, idx), ccol); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_fmadd_ps(a2, _mm256_permutevar8x32_ps(bcol, idx), ccol); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_fmadd_ps(a3, _mm256_permutevar8x32_ps(bcol, idx), ccol); idx = _mm256_add_epi32(idx, increament); a0 = _mm256_load_ps((float *)&left->vec[4][0]); a1 = _mm256_load_ps((float *)&left->vec[5][0]); a2 = _mm256_load_ps((float *)&left->vec[6][0]); a3 = _mm256_load_ps((float *)&left->vec[7][0]); ccol = _mm256_fmadd_ps(a0, _mm256_permutevar8x32_ps(bcol, idx), ccol); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_fmadd_ps(a1, _mm256_permutevar8x32_ps(bcol, idx), ccol); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_fmadd_ps(a2, _mm256_permutevar8x32_ps(bcol, idx), ccol); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_fmadd_ps(a3, _mm256_permutevar8x32_ps(bcol, idx), ccol); _mm256_store_ps((float *)&ans->vec[col][0], ccol); } } // ans = left * right + ans // right is not transposed DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, avx, fma) { static const __m256i idx0 = { 1, 1, 1, 1, 1, 1, 1, 1 }; for (int col = 0; col < 8; ++col) { register __m256 bcol = _mm256_load_ps((float *)&right->vec[col][0]); register __m256 ccol = _mm256_load_ps((float *)&ans->vec[col][0]); register __m256i idx = _mm256_setzero_si256(); register __m256i increament = _mm256_load_si256(&idx0); register __m256 a0 = _mm256_load_ps((float *)&left->vec[0][0]); register __m256 a1 = _mm256_load_ps((float *)&left->vec[1][0]); register __m256 a2 = _mm256_load_ps((float *)&left->vec[2][0]); register __m256 a3 = _mm256_load_ps((float *)&left->vec[3][0]); ccol = _mm256_add_ps(ccol, _mm256_mul_ps(a0, _mm256_permutevar8x32_ps(bcol, idx))); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_add_ps(ccol, _mm256_mul_ps(a1, _mm256_permutevar8x32_ps(bcol, idx))); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_add_ps(ccol, _mm256_mul_ps(a2, _mm256_permutevar8x32_ps(bcol, idx))); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_add_ps(ccol, _mm256_mul_ps(a3, _mm256_permutevar8x32_ps(bcol, idx))); idx = _mm256_add_epi32(idx, increament); a0 = _mm256_load_ps((float *)&left->vec[4][0]); a1 = _mm256_load_ps((float *)&left->vec[5][0]); a2 = _mm256_load_ps((float *)&left->vec[6][0]); a3 = _mm256_load_ps((float *)&left->vec[7][0]); ccol = _mm256_add_ps(ccol, _mm256_mul_ps(a0, _mm256_permutevar8x32_ps(bcol, idx))); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_add_ps(ccol, _mm256_mul_ps(a1, _mm256_permutevar8x32_ps(bcol, idx))); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_add_ps(ccol, _mm256_mul_ps(a2, _mm256_permutevar8x32_ps(bcol, idx))); idx = _mm256_add_epi32(idx, increament); ccol = _mm256_add_ps(ccol, _mm256_mul_ps(a3, _mm256_permutevar8x32_ps(bcol, idx))); _mm256_store_ps((float *)&ans->vec[col][0], ccol); } } DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, 512, fma) { const __m512i ones = _mm512_set1_epi32(1); for (int col = 0; col < 16; ++col) { register __m512 bcol = _mm512_load_ps((float *)&right->vec[col][0]); register __m512 ccol = _mm512_load_ps((float *)&ans->vec[col][0]); register __m512i idx = _mm512_setzero_si512(); register __m512i increament = _mm512_load_si512(&ones); register __m512 a0 = _mm512_load_ps((float *)&left->vec[0][0]); register __m512 a1 = _mm512_load_ps((float *)&left->vec[1][0]); register __m512 a2 = _mm512_load_ps((float *)&left->vec[2][0]); register __m512 a3 = _mm512_load_ps((float *)&left->vec[3][0]); register __m512 a4 = _mm512_load_ps((float *)&left->vec[4][0]); register __m512 a5 = _mm512_load_ps((float *)&left->vec[5][0]); register __m512 a6 = _mm512_load_ps((float *)&left->vec[6][0]); register __m512 a7 = _mm512_load_ps((float *)&left->vec[7][0]); ccol = _mm512_fmadd_ps(a0, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a1, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a2, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a3, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a4, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a5, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a6, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a7, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); a0 = _mm512_load_ps((float *)&left->vec[8][0]); a1 = _mm512_load_ps((float *)&left->vec[9][0]); a2 = _mm512_load_ps((float *)&left->vec[10][0]); a3 = _mm512_load_ps((float *)&left->vec[11][0]); a4 = _mm512_load_ps((float *)&left->vec[12][0]); a5 = _mm512_load_ps((float *)&left->vec[13][0]); a6 = _mm512_load_ps((float *)&left->vec[14][0]); a7 = _mm512_load_ps((float *)&left->vec[15][0]); ccol = _mm512_fmadd_ps(a0, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a1, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a2, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a3, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a4, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a5, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a6, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); ccol = _mm512_fmadd_ps(a7, _mm512_permutexvar_ps(idx, bcol), ccol); idx = _mm512_add_epi32(idx, increament); _mm512_store_ps((float *)&ans->vec[col][0], ccol); } } DECLARE_PNL_PROC_LEFT_RIGHT_ANS(double, sse, fma) { for (int col = 0; col < 4; ++col) { /* current answer column (two chunks: rows 0-1 and 2-3) */ register __m128d c_lo = _mm_castps_pd(PANEL_COL_VEC((*ans), col, 0)); register __m128d c_hi = _mm_castps_pd(PANEL_COL_VEC((*ans), col, 1)); /* right matrix column split into two chunks */ register __m128d b_lo = _mm_castps_pd(PANEL_COL_VEC((*right), col, 0)); register __m128d b_hi = _mm_castps_pd(PANEL_COL_VEC((*right), col, 1)); /* k = 0 */ { register __m128d a_lo = _mm_castps_pd(PANEL_COL_VEC((*left), 0, 0)); register __m128d a_hi = _mm_castps_pd(PANEL_COL_VEC((*left), 0, 1)); register __m128d b = _mm_shuffle_pd(b_lo, b_lo, 0x0); c_lo = _mm_add_pd(c_lo, _mm_mul_pd(a_lo, b)); c_hi = _mm_add_pd(c_hi, _mm_mul_pd(a_hi, b)); } /* k = 1 */ { register __m128d a_lo = _mm_castps_pd(PANEL_COL_VEC((*left), 1, 0)); register __m128d a_hi = _mm_castps_pd(PANEL_COL_VEC((*left), 1, 1)); register __m128d b = _mm_shuffle_pd(b_lo, b_lo, 0x3); c_lo = _mm_add_pd(c_lo, _mm_mul_pd(a_lo, b)); c_hi = _mm_add_pd(c_hi, _mm_mul_pd(a_hi, b)); } /* k = 2 */ { register __m128d a_lo = _mm_castps_pd(PANEL_COL_VEC((*left), 2, 0)); register __m128d a_hi = _mm_castps_pd(PANEL_COL_VEC((*left), 2, 1)); register __m128d b = _mm_shuffle_pd(b_hi, b_hi, 0x0); c_lo = _mm_add_pd(c_lo, _mm_mul_pd(a_lo, b)); c_hi = _mm_add_pd(c_hi, _mm_mul_pd(a_hi, b)); } /* k = 3 */ { register __m128d a_lo = _mm_castps_pd(PANEL_COL_VEC((*left), 3, 0)); register __m128d a_hi = _mm_castps_pd(PANEL_COL_VEC((*left), 3, 1)); register __m128d b = _mm_shuffle_pd(b_hi, b_hi, 0x3); c_lo = _mm_add_pd(c_lo, _mm_mul_pd(a_lo, b)); c_hi = _mm_add_pd(c_hi, _mm_mul_pd(a_hi, b)); } PANEL_COL_VEC((*ans), col, 0) = _mm_castpd_ps(c_lo); PANEL_COL_VEC((*ans), col, 1) = _mm_castpd_ps(c_hi); } } DECLARE_PNL_PROC_LEFT_RIGHT_ANS(double, avx, fma) { register __m256d a0 = _mm256_load_pd((double *)&left->vec[0][0]); register __m256d a1 = _mm256_load_pd((double *)&left->vec[1][0]); register __m256d a2 = _mm256_load_pd((double *)&left->vec[2][0]); register __m256d a3 = _mm256_load_pd((double *)&left->vec[3][0]); register __m256d c0 = _mm256_load_pd((double *)&ans->vec[0][0]); register __m256d c1 = _mm256_load_pd((double *)&ans->vec[1][0]); register __m256d c2 = _mm256_load_pd((double *)&ans->vec[2][0]); register __m256d c3 = _mm256_load_pd((double *)&ans->vec[3][0]); register __m256d bcol = _mm256_load_pd((double *)&right->vec[0][0]); c0 = _mm256_add_pd(c0, _mm256_mul_pd(a0, _mm256_permute_pd(bcol, 0x0))); c0 = _mm256_add_pd(c0, _mm256_mul_pd(a1, _mm256_permute_pd(bcol, 0x3))); c0 = _mm256_add_pd(c0, _mm256_mul_pd(a2, _mm256_permute_pd(bcol, 0x0))); c0 = _mm256_add_pd(c0, _mm256_mul_pd(a3, _mm256_permute_pd(bcol, 0x3))); _mm256_store_pd((double *)&ans->vec[0][0], c0); bcol = _mm256_load_pd((double *)&right->vec[0][0]); c1 = _mm256_add_pd(c1, _mm256_mul_pd(a0, _mm256_permute_pd(bcol, 0x0))); c1 = _mm256_add_pd(c1, _mm256_mul_pd(a1, _mm256_permute_pd(bcol, 0x3))); c1 = _mm256_add_pd(c1, _mm256_mul_pd(a2, _mm256_permute_pd(bcol, 0x0))); c1 = _mm256_add_pd(c1, _mm256_mul_pd(a3, _mm256_permute_pd(bcol, 0x3))); _mm256_store_pd((double *)&ans->vec[1][0], c1); bcol = _mm256_load_pd((double *)&right->vec[0][0]); c2 = _mm256_add_pd(c2, _mm256_mul_pd(a0, _mm256_permute_pd(bcol, 0x0))); c2 = _mm256_add_pd(c2, _mm256_mul_pd(a1, _mm256_permute_pd(bcol, 0x3))); c2 = _mm256_add_pd(c2, _mm256_mul_pd(a2, _mm256_permute_pd(bcol, 0x0))); c2 = _mm256_add_pd(c2, _mm256_mul_pd(a3, _mm256_permute_pd(bcol, 0x3))); _mm256_store_pd((double *)&ans->vec[2][0], c2); bcol = _mm256_load_pd((double *)&right->vec[0][0]); c3 = _mm256_add_pd(c3, _mm256_mul_pd(a0, _mm256_permute_pd(bcol, 0x0))); c3 = _mm256_add_pd(c3, _mm256_mul_pd(a1, _mm256_permute_pd(bcol, 0x3))); c3 = _mm256_add_pd(c3, _mm256_mul_pd(a2, _mm256_permute_pd(bcol, 0x0))); c3 = _mm256_add_pd(c3, _mm256_mul_pd(a3, _mm256_permute_pd(bcol, 0x3))); _mm256_store_pd((double *)&ans->vec[3][0], c3); } DECLARE_PNL_PROC_LEFT_RIGHT_ANS(double, avx_fma, fma) { register __m256d a0 = _mm256_load_pd((double *)&left->vec[0][0]); register __m256d a1 = _mm256_load_pd((double *)&left->vec[1][0]); register __m256d a2 = _mm256_load_pd((double *)&left->vec[2][0]); register __m256d a3 = _mm256_load_pd((double *)&left->vec[3][0]); register __m256d c0 = _mm256_load_pd((double *)&ans->vec[0][0]); register __m256d c1 = _mm256_load_pd((double *)&ans->vec[1][0]); register __m256d c2 = _mm256_load_pd((double *)&ans->vec[2][0]); register __m256d c3 = _mm256_load_pd((double *)&ans->vec[3][0]); register __m256d bcol = _mm256_load_pd((double *)&right->vec[0][0]); c0 = _mm256_fmadd_pd(a0, _mm256_permute_pd(bcol, 0x0), c0); c0 = _mm256_fmadd_pd(a1, _mm256_permute_pd(bcol, 0x3), c0); c0 = _mm256_fmadd_pd(a2, _mm256_permute_pd(bcol, 0x0), c0); c0 = _mm256_fmadd_pd(a3, _mm256_permute_pd(bcol, 0x3), c0); _mm256_store_pd((double *)&ans->vec[0][0], c0); bcol = _mm256_load_pd((double *)&right->vec[0][0]); c1 = _mm256_fmadd_pd(a0, _mm256_permute_pd(bcol, 0x0), c1); c1 = _mm256_fmadd_pd(a1, _mm256_permute_pd(bcol, 0x3), c1); c1 = _mm256_fmadd_pd(a2, _mm256_permute_pd(bcol, 0x0), c1); c1 = _mm256_fmadd_pd(a3, _mm256_permute_pd(bcol, 0x3), c1); _mm256_store_pd((double *)&ans->vec[1][0], c1); bcol = _mm256_load_pd((double *)&right->vec[0][0]); c2 = _mm256_fmadd_pd(a0, _mm256_permute_pd(bcol, 0x0), c2); c2 = _mm256_fmadd_pd(a1, _mm256_permute_pd(bcol, 0x3), c2); c2 = _mm256_fmadd_pd(a2, _mm256_permute_pd(bcol, 0x0), c2); c2 = _mm256_fmadd_pd(a3, _mm256_permute_pd(bcol, 0x3), c2); _mm256_store_pd((double *)&ans->vec[2][0], c2); bcol = _mm256_load_pd((double *)&right->vec[0][0]); c3 = _mm256_fmadd_pd(a0, _mm256_permute_pd(bcol, 0x0), c3); c3 = _mm256_fmadd_pd(a1, _mm256_permute_pd(bcol, 0x3), c3); c3 = _mm256_fmadd_pd(a2, _mm256_permute_pd(bcol, 0x0), c3); c3 = _mm256_fmadd_pd(a3, _mm256_permute_pd(bcol, 0x3), c3); _mm256_store_pd((double *)&ans->vec[3][0], c3); } // ans = left * right + ans // right is not transposed DECLARE_PNL_PROC_LEFT_RIGHT_ANS(double, 512, fma) { const __m512i ones = _mm512_set1_epi32(1); register __m512d a0 = _mm512_load_pd((double *)&left->vec[0][0]); register __m512d a1 = _mm512_load_pd((double *)&left->vec[1][0]); register __m512d a2 = _mm512_load_pd((double *)&left->vec[2][0]); register __m512d a3 = _mm512_load_pd((double *)&left->vec[3][0]); register __m512d a4 = _mm512_load_pd((double *)&left->vec[4][0]); register __m512d a5 = _mm512_load_pd((double *)&left->vec[5][0]); register __m512d a6 = _mm512_load_pd((double *)&left->vec[6][0]); register __m512d a7 = _mm512_load_pd((double *)&left->vec[7][0]); for (int col = 0; col < 8; ++col) { register __m512d bcol = _mm512_load_pd((double *)&right->vec[col][0]); register __m512d ccol = _mm512_load_pd((double *)&ans->vec[col][0]); register __m512i idx = _mm512_setzero_si512(); register __m512i increament = _mm512_load_si512(&ones); ccol = _mm512_fmadd_pd(a0, _mm512_permutexvar_pd(idx, bcol), ccol); idx = _mm512_add_epi64(idx, increament); ccol = _mm512_fmadd_pd(a1, _mm512_permutexvar_pd(idx, bcol), ccol); idx = _mm512_add_epi64(idx, increament); ccol = _mm512_fmadd_pd(a2, _mm512_permutexvar_pd(idx, bcol), ccol); idx = _mm512_add_epi64(idx, increament); ccol = _mm512_fmadd_pd(a3, _mm512_permutexvar_pd(idx, bcol), ccol); idx = _mm512_add_epi64(idx, increament); ccol = _mm512_fmadd_pd(a4, _mm512_permutexvar_pd(idx, bcol), ccol); idx = _mm512_add_epi64(idx, increament); ccol = _mm512_fmadd_pd(a5, _mm512_permutexvar_pd(idx, bcol), ccol); idx = _mm512_add_epi64(idx, increament); ccol = _mm512_fmadd_pd(a6, _mm512_permutexvar_pd(idx, bcol), ccol); idx = _mm512_add_epi64(idx, increament); ccol = _mm512_fmadd_pd(a7, _mm512_permutexvar_pd(idx, bcol), ccol); idx = _mm512_add_epi64(idx, increament); _mm512_store_pd((double *)&ans->vec[col][0], ccol); } }