I in the process of making some 2D/planar image data resampler and current highest performance approach of 2D convolution need to perform shift of long float32 vector between AVX (ymm or zmm) registers at fixed number of float32 (4 byte) values. The shift distance is equal to multiplication ratio of image enlargment, so for integer ratio I need to shift to 2,3,4,5,6 float32 in between long number of AVX registers. The shifting on 8 and 16 is just a renaming of registers.
The fastest current code for shift to 4 float32 uses rare AVX instruction
_mm256_permute2f128_ps. It allow to shift and exchange data between 2 ymm registers.
So current shift to 4 float32 in between 7 ymm registers looks like
_mm_store_ps(pfProc, _mm256_castps256_ps128(my_ymm0));
my_ymm0 = _mm256_permute2f128_ps(my_ymm0, my_ymm1, 33);
my_ymm1 = _mm256_permute2f128_ps(my_ymm1, my_ymm2, 33);
my_ymm2 = _mm256_permute2f128_ps(my_ymm2, my_ymm3, 33);
my_ymm3 = _mm256_permute2f128_ps(my_ymm3, my_ymm4, 33);
my_ymm4 = _mm256_permute2f128_ps(my_ymm4, my_ymm5, 33);
my_ymm5 = _mm256_permute2f128_ps(my_ymm5, my_ymm6, 33);
my_ymm6 = _mm256_permute2f128_ps(my_ymm6, my_ymm6, 49);
my_ymm6 = _mm256_insertf128_ps(my_ymm6, *(__m128*)(pfProc + 56), 1);
with storing out 4 float result at first, shift and load 4 float to the end. With a large number of fma instructions this multi-sample convolution approach reach about 30..35% of theoretical FMA performance of 4 core intel 9th generation CPU with 4 FMA512 units.
Unfortunately I can not found 1 instruction to make shift to other number of floats in the instruction sets before AVX512VL.
Currently found approach on before_AVX512VL for shifting to 2 floats uses permute for shift inside ymm and blend to transfer shifted-out number of floats to another ymm:
const register __m256i my_ymm8_main_circ = _mm256_set_epi32(1, 0, 7, 6, 5, 4, 3, 2); // main circulating const
my_ymm2 = _mm256_permutevar8x32_ps(my_ymm2, my_ymm8_main_circ); // circulate by 2 ps to the left
my_ymm3 = _mm256_permutevar8x32_ps(my_ymm3, my_ymm8_main_circ); // circulate by 2 ps to the left
my_ymm2 = _mm256_blend_ps(my_ymm2, my_ymm3, 192); // copy higher 2 floats
my_ymm4 = _mm256_permutevar8x32_ps(my_ymm4, my_ymm8_main_circ); // circulate by 2 ps to the left
my_ymm3 = _mm256_blend_ps(my_ymm3, my_ymm4, 192); // copy higher 2 floats
my_ymm5 = _mm256_permutevar8x32_ps(my_ymm5, my_ymm8_main_circ); // circulate by 2 ps to the left
my_ymm4 = _mm256_blend_ps(my_ymm4, my_ymm5, 192); // copy higher 2 floats
my_ymm6 = _mm256_permutevar8x32_ps(my_ymm6, my_ymm8_main_circ); // circulate by 2 ps to the left
my_ymm5 = _mm256_blend_ps(my_ymm5, my_ymm6, 192); // copy higher 2 floats
my_ymm7 = _mm256_permutevar8x32_ps(my_ymm7, my_ymm8_main_circ); // circulate by 2 ps to the left
my_ymm6 = _mm256_blend_ps(my_ymm6, my_ymm7, 192); // copy higher 2 floats
It works visibly slower and uses very large immediate (sized as large as ymm 256 bit register) to set only a few control bits for permutation. That do not allow do use more ymm regs for data storing and cause compiler to load/store immediate and slow down execution.
The other approach is to store (to L1d cache) the long vector from 5..7 ymm regs and load back with shifted address. But it looks even slower because of large buck memory read-write operation and I can not be sure if it will not be transfered via all cache hierarchy causing memory controller to slow down the main memory load/store streams in multi core CPU.
With AVX512VL instruction set I see the 2-sources permutation instruction
_mm512_permutex2var_ps (also uses large zmm immediate to select float32 for permutation).
It can replace _mm256_permute2f128_ps for any number of float32 taken from 1 and from 2nd zmm registers so allow to shift to any number of float32 between registers in 1 instruction.
But AVX512VL CPUs are not in common use at endusers in even 2021. May be in 10 years later will. But not now.
May be other approaches exist for non-AVX512VL capable CPUs for faster shift of float32 long vectors beween a number of ymm (and may be zmm) registers at any step measured in ps floats ?
I got several questions about the whole implemented algorithm using this long vector shift. One of the most short 2D/planar upsizer to 4x function listing is:
void JincResize::KernelRowAll_avx2_mul4_taps4_cb(unsigned char* src, int iSrcStride, int iInpWidth, int iInpHeight, unsigned char* dst, int iDstStride, const int ciMul, const int ciTaps, const int ciKS)
{
// current input plane sizes
iWidthEl = iInpWidth + 64/*2 * iKernelSize*/;
iHeightEl = iInpHeight + 64/*2 * iKernelSize*/;
const int k_col8 = 32/*iKernelSize*/ - (32/*iKernelSize*/ % 8);
const int col8 = iWidthEl - 4/*iTaps*/ - ((iWidthEl - 4/*iTaps*/) % 8);
float* pfCurrKernel = g_pfKernel;
int64_t col;
memset(pfFilteredCirculatingBuf, 0, iWidthEl * 32 * 4/*iKernelSize * iMul*/ * sizeof(float));
// no MT version
for (int64_t row = 4/*iTaps*/; row < iHeightEl - 4/*iTaps*/; row++) // input lines counter
{
// start all row-only dependent ptrs here
float* pfInpRowSamplesFloatBufStart = pfInpFloatRow;
GetInpElRowAsFloat_avx2(row, iInpHeight, iInpWidth, src, iSrcStride, pfInpRowSamplesFloatBufStart, 32/*iKernelSize*/);
for (col = 4/*iTaps*/; col < col8; col += 8) // input cols counter
{
float* pfColStart = pfInpRowSamplesFloatBufStart + col;
float* pfCurrKernel_pos = pfCurrKernel;
for (int64_t k_row = 0; k_row < iKernelSize; k_row++)
{
float* pfProc = vpfRowsPointers[k_row] + col * 4/*iMul*/;
__m256 my_ymm0, my_ymm1, my_ymm2, my_ymm3, my_ymm4, my_ymm5, my_ymm6, my_ymm7; // out samples
__m256 my_ymm8, my_ymm9, my_ymm10, my_ymm11; // inp samples
__m256 my_ymm12, my_ymm13, my_ymm14, my_ymm15; // kernel samples
my_ymm12 = _mm256_load_ps(pfCurrKernel_pos);
my_ymm13 = _mm256_load_ps(pfCurrKernel_pos + 8);
my_ymm14 = _mm256_load_ps(pfCurrKernel_pos + 16);
my_ymm15 = _mm256_load_ps(pfCurrKernel_pos + 24);
my_ymm0 = _mm256_load_ps(pfProc);
my_ymm1 = _mm256_load_ps(pfProc + 8);
my_ymm2 = _mm256_load_ps(pfProc + 16);
my_ymm3 = _mm256_load_ps(pfProc + 24);
my_ymm4 = _mm256_load_ps(pfProc + 32);
my_ymm5 = _mm256_load_ps(pfProc + 40);
my_ymm6 = _mm256_load_ps(pfProc + 48);
my_ymm8 = _mm256_broadcast_ss(pfColStart + 0); // 1
my_ymm9 = _mm256_broadcast_ss(pfColStart + 2); // 3
my_ymm10 = _mm256_broadcast_ss(pfColStart + 4); // 5
my_ymm11 = _mm256_broadcast_ss(pfColStart + 6); // 7
// 1st sample
my_ymm0 = _mm256_fmadd_ps(my_ymm12, my_ymm8, my_ymm0);
my_ymm1 = _mm256_fmadd_ps(my_ymm13, my_ymm8, my_ymm1);
my_ymm2 = _mm256_fmadd_ps(my_ymm14, my_ymm8, my_ymm2);
my_ymm3 = _mm256_fmadd_ps(my_ymm15, my_ymm8, my_ymm3);
// 3rd sample
my_ymm1 = _mm256_fmadd_ps(my_ymm12, my_ymm9, my_ymm1);
my_ymm2 = _mm256_fmadd_ps(my_ymm13, my_ymm9, my_ymm2);
my_ymm3 = _mm256_fmadd_ps(my_ymm14, my_ymm9, my_ymm3);
my_ymm4 = _mm256_fmadd_ps(my_ymm15, my_ymm9, my_ymm4);
// 5th sample
my_ymm2 = _mm256_fmadd_ps(my_ymm12, my_ymm10, my_ymm2);
my_ymm3 = _mm256_fmadd_ps(my_ymm13, my_ymm10, my_ymm3);
my_ymm4 = _mm256_fmadd_ps(my_ymm14, my_ymm10, my_ymm4);
my_ymm5 = _mm256_fmadd_ps(my_ymm15, my_ymm10, my_ymm5);
// 7th sample
my_ymm3 = _mm256_fmadd_ps(my_ymm12, my_ymm11, my_ymm3);
my_ymm4 = _mm256_fmadd_ps(my_ymm13, my_ymm11, my_ymm4);
my_ymm5 = _mm256_fmadd_ps(my_ymm14, my_ymm11, my_ymm5);
my_ymm6 = _mm256_fmadd_ps(my_ymm15, my_ymm11, my_ymm6);
_mm_store_ps(pfProc, _mm256_castps256_ps128(my_ymm0));
my_ymm0 = _mm256_permute2f128_ps(my_ymm0, my_ymm1, 33);
my_ymm1 = _mm256_permute2f128_ps(my_ymm1, my_ymm2, 33);
my_ymm2 = _mm256_permute2f128_ps(my_ymm2, my_ymm3, 33);
my_ymm3 = _mm256_permute2f128_ps(my_ymm3, my_ymm4, 33);
my_ymm4 = _mm256_permute2f128_ps(my_ymm4, my_ymm5, 33);
my_ymm5 = _mm256_permute2f128_ps(my_ymm5, my_ymm6, 33);
my_ymm6 = _mm256_permute2f128_ps(my_ymm6, my_ymm6, 49);
my_ymm6 = _mm256_insertf128_ps(my_ymm6, *(__m128*)(pfProc + 56), 1);
// even samples
my_ymm8 = _mm256_broadcast_ss(pfColStart + 1); // 2
my_ymm9 = _mm256_broadcast_ss(pfColStart + 3); // 4
my_ymm10 = _mm256_broadcast_ss(pfColStart + 5); // 6
my_ymm11 = _mm256_broadcast_ss(pfColStart + 7); // 8
// 2nd sample
my_ymm0 = _mm256_fmadd_ps(my_ymm12, my_ymm8, my_ymm0);
my_ymm1 = _mm256_fmadd_ps(my_ymm13, my_ymm8, my_ymm1);
my_ymm2 = _mm256_fmadd_ps(my_ymm14, my_ymm8, my_ymm2);
my_ymm3 = _mm256_fmadd_ps(my_ymm15, my_ymm8, my_ymm3);
// 4th sample
my_ymm1 = _mm256_fmadd_ps(my_ymm12, my_ymm9, my_ymm1);
my_ymm2 = _mm256_fmadd_ps(my_ymm13, my_ymm9, my_ymm2);
my_ymm3 = _mm256_fmadd_ps(my_ymm14, my_ymm9, my_ymm3);
my_ymm4 = _mm256_fmadd_ps(my_ymm15, my_ymm9, my_ymm4);
// 6th sample
my_ymm2 = _mm256_fmadd_ps(my_ymm12, my_ymm10, my_ymm2);
my_ymm3 = _mm256_fmadd_ps(my_ymm13, my_ymm10, my_ymm3);
my_ymm4 = _mm256_fmadd_ps(my_ymm14, my_ymm10, my_ymm4);
my_ymm5 = _mm256_fmadd_ps(my_ymm15, my_ymm10, my_ymm5);
// 7th sample
my_ymm3 = _mm256_fmadd_ps(my_ymm12, my_ymm11, my_ymm3);
my_ymm4 = _mm256_fmadd_ps(my_ymm13, my_ymm11, my_ymm4);
my_ymm5 = _mm256_fmadd_ps(my_ymm14, my_ymm11, my_ymm5);
my_ymm6 = _mm256_fmadd_ps(my_ymm15, my_ymm11, my_ymm6);
_mm256_store_ps(pfProc + 4, my_ymm0);
_mm256_store_ps(pfProc + 12, my_ymm1);
_mm256_store_ps(pfProc + 20, my_ymm2);
_mm256_store_ps(pfProc + 28, my_ymm3);
_mm256_store_ps(pfProc + 36, my_ymm4);
_mm256_store_ps(pfProc + 44, my_ymm5);
_mm256_store_ps(pfProc + 52, my_ymm6);
pfCurrKernel_pos += 32/*iKernelSize*/; // point to next kernel row now
} // k_row
} // col
// need to process last up to 7 cols separately...
for (col = col8 + 4/*iTaps*/; col < iWidthEl - 4/*iTaps*/; col++) // input cols counter
{
float* pfCurrKernel_pos = pfCurrKernel;
float* pfProc;
for (int64_t k_row = 0; k_row < 32/*iKernelSize*/; k_row++)
{
pfProc = vpfRowsPointers[k_row] + col * 4/*iMul*/;
for (int64_t k_col = 0; k_col < k_col8; k_col += 8)
{
*(__m256*)(pfProc + k_col) = _mm256_fmadd_ps(*(__m256*)(pfCurrKernel_pos + k_col), _mm256_broadcast_ss(&pfInpRowSamplesFloatBufStart[col]), *(__m256*)(pfProc + k_col));
}
// need to process last (up to 7) floats separately..
for (int64_t k_col = k_col8; k_col < 32/*iKernelSize*/; ++k_col)
{
pfProc[k_col] += (pfCurrKernel_pos[k_col] * pfInpRowSamplesFloatBufStart[col]);
}
pfCurrKernel_pos += 32/*iKernelSize*/; // point to next kernel row now
} // k_row
} // col
int iOutStartRow = (row - (4 + 32/*iTaps + iKernelSize*/)) * 4/*iMul*/;
//iMul rows ready - output result, skip iKernelSize+iTaps rows from beginning
if (iOutStartRow >= 0 && iOutStartRow < (iInpHeight) * 4/*iMul*/)
{
ConvertiMulRowsToInt_avx2(vpfRowsPointers, iInpWidth, iOutStartRow, dst, iDstStride);
}
// circulate pointers to iMul rows upper
std::rotate(vpfRowsPointers.begin(), vpfRowsPointers.begin() + 4/*iMul*/, vpfRowsPointers.end());
// clear last iMul rows
for (int i = 32 - 4/*iKernelSize - iMul*/; i < 32/*iKernelSize*/; i++)
{
memset(vpfRowsPointers[i], 0, iWidthEl * 4/*iMul*/ * sizeof(float));
}
} // row
}
As asked adding completely non-vectorized C-version of resampler:
void JincResize::KernelRowAll_c_mul_cb(unsigned char *src, int iSrcStride, int iInpWidth, int iInpHeight, unsigned char *dst, int iDstStride)
{
iWidthEl = iInpWidth + 2 * iKernelSize;
iHeightEl = iInpHeight + 2 * iKernelSize;
float* pfCurrKernel = g_pfKernel;
memset(pfFilteredCirculatingBuf, 0, iWidthEl * iKernelSize * iMul * num_threads * sizeof(float));
// single threaded for now
for (int64_t row = iTaps; row < iHeightEl - iTaps; row++) // input lines counter
{
float* pfInpRowSamplesFloatBufStart = pfInpFloatRow;
(this->*GetInpElRowAsFloat)(row, iInpHeight, iInpWidth, src, iSrcStride, pfInpRowSamplesFloatBufStart);
for (int64_t col = iTaps; col < iWidthEl - iTaps; col++) // input cols counter
{
float fInpSample = pfInpRowSamplesFloatBufStart[col];
float* pfCurrKernel_pos = pfCurrKernel;
float* pfProc;
for (int64_t k_row = 0; k_row < iKernelSize; k_row++)
{
pfProc = vpfRowsPointers[k_row] + col * iMul;
for (int64_t k_col = 0; k_col < iKernelSize; k_col++)
{
pfProc[k_col] += pfCurrKernel_pos[k_col] * fInpSample;
} // k_col
pfCurrKernel_pos += iKernelSize; // point to next kernel row now
} // k_row
} // col
int iOutStartRow = (row - (iTaps + iKernelSize))*iMul;
//iMul rows ready - output result, skip iKernelSize+iTaps rows from beginning
if (iOutStartRow >= 0 && iOutStartRow < (iInpHeight)*iMul)
{
ConvertNRowsToInt_c(vpfRowsPointers, iInpWidth, iOutStartRow, dst, iDstStride, iMul);
}
// circulate pointers to iMul rows upper
std::rotate(vpfRowsPointers.begin(), vpfRowsPointers.begin() + iMul, vpfRowsPointers.end());
// clear last iMul rows
for (int i = iKernelSize - iMul; i < iKernelSize; i++)
{
memset(vpfRowsPointers[i], 0, iWidthEl*iMul * sizeof(float));
}
}
}