Floating point numbers and unsigned bytes are both used in graphics. Generally 0..1 for float and obviously 0..255 for unsigned byte.
It's nice when we can stick to one or the other, but sometimes that doesn't work and we have to convert between them.
There are details. Once again there are always details. And I've seen people get them wrong. Often I've seen this implemented with truncation rather than rounding - which makes me sad. I've also seen a direct cast to byte, without limiting to zero and 255 - which makes me sad.
Pretty much all of the float SSE instructions are affected by the rounding mode. Pretty much nobody changes these bits without changing them back when they're done. Setting the bits is a little expensive - pretty much everything has to stop in order to do this. But amortized over a function call it's not bad.
__forceinline unsigned char Saturate_Int32_To_UnsignedInt8(const int x)
{
return std::min(255, std::max(0, x));
}
__forceinline unsigned char Saturate_float_To_UnsignedInt8(const float x)
{
// since we can't set the rounding mode for simple ANSI, we'll do the normal thing and add 0.5f and truncate towards zero.
int ic = (int)((255 * x)+0.5f);
return Saturate_Int32_To_UnsignedInt8(ic);
}
void convert_floats_to_bytes(const float *pF, unsigned char *pUC, const size_t cF)
{
// setting the rounding mode is generally unnecessary... But to be safe....
const unsigned int uiStoredRoundingMode = _MM_GET_ROUNDING_MODE();
_MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST);
size_t cFRemaining = cF;
while (cFRemaining & 3)
{
cFRemaining--;
pUC[cFRemaining] = Saturate_float_To_UnsignedInt8(pF[cFRemaining]);
}
if (cFRemaining > 0)
{
__m128 mScale = _mm_set1_ps(255);
if (cFRemaining & 4)
{
cFRemaining -= 4;
__m128 mF = _mm_loadu_ps(&pF[cFRemaining]);
__m128 mFScaled = _mm_mul_ps(mF, mScale);
__m128i mIScaled_epi32 = _mm_cvtps_epi32(mFScaled);
__m128i mIScaled_epi16 = _mm_packs_epi32(mIScaled_epi32, _mm_setzero_si128());
__m128i mIScaled_epu8 = _mm_packus_epi16(mIScaled_epi16, _mm_setzero_si128());
int iConvertedOut;
iConvertedOut = _mm_cvtsi128_si32(mIScaled_epu8);
*((int *)&pUC[cFRemaining]) = iConvertedOut;
}
while (cFRemaining > 0)
{
cFRemaining -= 4;
__m128 mFHi = _mm_loadu_ps(&pF[cFRemaining]);
cFRemaining -= 4;
__m128 mFLo = _mm_loadu_ps(&pF[cFRemaining]);
__m128 mFScaledHi = _mm_mul_ps(mFHi, mScale);
__m128 mFScaledLo = _mm_mul_ps(mFLo, mScale);
__m128i mIScaledHi_epi32 = _mm_cvtps_epi32(mFScaledHi);
__m128i mIScaledLo_epi32 = _mm_cvtps_epi32(mFScaledLo);
__m128i mIScaled_epi16 = _mm_packs_epi32(mIScaledLo_epi32, mIScaledHi_epi32);
__m128i mIScaled_epu8 = _mm_packus_epi16(mIScaled_epi16, _mm_setzero_si128());
__int64 iConvertedOut;
#if defined (_M_X64)
iConvertedOut = _mm_cvtsi128_si64(mIScaled_epu8);
#else
iConvertedOut = mIScaled_epu8.m128i_i64[0];
#endif
*((__int64 *)&pUC[cFRemaining]) = iConvertedOut;
}
}
_MM_SET_ROUNDING_MODE(uiStoredRoundingMode);
}
This doesn't make me sad.
Thursday, March 23, 2017
Tuesday, March 21, 2017
a bitonic sorting network for 16 signed short in SSE2
We saw a sorting network for 8 float numbers. These values fit nicely into two 128bit registers.
Twice as many signed shorts can be fit into the same space. And SSE2 has intrinsics for maximum and minimum for 8 signed short numbers at a time using _mm_min_epi16 and _mm_max_epi16.
While the bitonic sorting network for 16 values requires 10 min/max pairs rather than the 9 required in the network shown in Knuth, the shuffles are much easier for the plain vanilla bitonic sorting network.
A lot of shuffling has to happen to make this all work. The unpack operations will happen in pairs:
void unpack_epi16_ab(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
__m128i aa = _mm_unpacklo_epi16(a, b);
__m128i bb = _mm_unpackhi_epi16(a, b);
aOut = aa;
bOut = bb;
}
void unpack_epi32_ab(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
__m128i aa = _mm_unpacklo_epi32(a, b);
__m128i bb = _mm_unpackhi_epi32(a, b);
aOut = aa;
bOut = bb;
}
void unpack_epi64_ab(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
__m128i aa = _mm_unpacklo_epi64(a, b);
__m128i bb = _mm_unpackhi_epi64(a, b);
aOut = aa;
bOut = bb;
}
Note that it is safe to call all of these with the same variables used as input and output.
void BitonicSortingNetwork_epi16_SSE2(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
__m128i a1 = a;
__m128i b1 = b;
__m128i a1Min = _mm_min_epi16(a1, b1);
__m128i b1Max = _mm_max_epi16(a1, b1);
__m128i a2 = a1Min;
__m128i b2 = _mm_shufflelo_epi16(b1Max, _MM_SHUFFLE(2, 3, 0, 1));
b2 = _mm_shufflehi_epi16(b2, _MM_SHUFFLE(2, 3, 0, 1));
__m128i a2Min = _mm_min_epi16(a2, b2);
__m128i b2Max = _mm_max_epi16(a2, b2);
__m128i a3, b3;
unpack_epi16_ab(a2Min, b2Max, a3, b3);
unpack_epi16_ab(a3, b3, a3, b3);
unpack_epi16_ab(a3, b3, a3, b3);
__m128i a3Min = _mm_min_epi16(a3, b3);
__m128i b3Max = _mm_max_epi16(a3, b3);
__m128i a4, b4;
unpack_epi32_ab(a3Min, b3Max, a4, b4);
a4 = _mm_shufflelo_epi16(a4, _MM_SHUFFLE(0, 1, 2, 3));
b4 = _mm_shufflehi_epi16(b4, _MM_SHUFFLE(0, 1, 2, 3));
__m128i a4Min = _mm_min_epi16(a4, b4);
__m128i b4Max = _mm_max_epi16(a4, b4);
__m128i a5, b5;
unpack_epi16_ab(a4Min, b4Max, a5, b5);
unpack_epi64_ab(a5, b5, a5, b5);
a5 = _mm_shuffle_epi32(a5, _MM_SHUFFLE(2, 3, 0, 1));
__m128i a5Min = _mm_min_epi16(a5, b5);
__m128i b5Max = _mm_max_epi16(a5, b5);
__m128i a6, b6;
unpack_epi16_ab(a5Min, b5Max, a6, b6);
unpack_epi16_ab(a6, b6, a6, b6);
__m128i a6Min = _mm_min_epi16(a6, b6);
__m128i b6Max = _mm_max_epi16(a6, b6);
__m128i a7, b7;
b7 = _mm_shuffle_epi32(b6Max, _MM_SHUFFLE(1, 0, 3, 2));
b7 = _mm_shufflehi_epi16(b7, _MM_SHUFFLE(0, 1, 2, 3));
a7 = _mm_shufflelo_epi16(a6Min, _MM_SHUFFLE(0, 1, 2, 3));
__m128i a7Min = _mm_min_epi16(a7, b7);
__m128i b7Max = _mm_max_epi16(a7, b7);
__m128i a8, b8;
unpack_epi16_ab(a7Min, b7Max, b8, a8);
a8 = _mm_shuffle_epi32(a8, _MM_SHUFFLE(0, 1, 2, 3));
__m128i a8Min = _mm_min_epi16(a8, b8);
__m128i b8Max = _mm_max_epi16(a8, b8);
__m128i a9, b9;
unpack_epi16_ab(a8Min, b8Max, a9, b9);
__m128i a9Min = _mm_min_epi16(a9, b9);
__m128i b9Max = _mm_max_epi16(a9, b9);
__m128i a10, b10;
unpack_epi16_ab(a9Min, b9Max, a10, b10);
__m128i a10Min = _mm_min_epi16(a10, b10);
__m128i b10Max = _mm_max_epi16(a10, b10);
unpack_epi16_ab(a10Min, b10Max, aOut, bOut);
}
This is verified using the standard 0/1 technique described in Knuth.
The compiled code uses only five registers, so avoids register spills.
I am still not sure if this was worth the effort. But it's done.
Twice as many signed shorts can be fit into the same space. And SSE2 has intrinsics for maximum and minimum for 8 signed short numbers at a time using _mm_min_epi16 and _mm_max_epi16.
While the bitonic sorting network for 16 values requires 10 min/max pairs rather than the 9 required in the network shown in Knuth, the shuffles are much easier for the plain vanilla bitonic sorting network.
A lot of shuffling has to happen to make this all work. The unpack operations will happen in pairs:
void unpack_epi16_ab(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
__m128i aa = _mm_unpacklo_epi16(a, b);
__m128i bb = _mm_unpackhi_epi16(a, b);
aOut = aa;
bOut = bb;
}
void unpack_epi32_ab(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
__m128i aa = _mm_unpacklo_epi32(a, b);
__m128i bb = _mm_unpackhi_epi32(a, b);
aOut = aa;
bOut = bb;
}
void unpack_epi64_ab(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
__m128i aa = _mm_unpacklo_epi64(a, b);
__m128i bb = _mm_unpackhi_epi64(a, b);
aOut = aa;
bOut = bb;
}
Note that it is safe to call all of these with the same variables used as input and output.
void BitonicSortingNetwork_epi16_SSE2(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
__m128i a1 = a;
__m128i b1 = b;
__m128i a1Min = _mm_min_epi16(a1, b1);
__m128i b1Max = _mm_max_epi16(a1, b1);
__m128i a2 = a1Min;
__m128i b2 = _mm_shufflelo_epi16(b1Max, _MM_SHUFFLE(2, 3, 0, 1));
b2 = _mm_shufflehi_epi16(b2, _MM_SHUFFLE(2, 3, 0, 1));
__m128i a2Min = _mm_min_epi16(a2, b2);
__m128i b2Max = _mm_max_epi16(a2, b2);
__m128i a3, b3;
unpack_epi16_ab(a2Min, b2Max, a3, b3);
unpack_epi16_ab(a3, b3, a3, b3);
unpack_epi16_ab(a3, b3, a3, b3);
__m128i a3Min = _mm_min_epi16(a3, b3);
__m128i b3Max = _mm_max_epi16(a3, b3);
__m128i a4, b4;
unpack_epi32_ab(a3Min, b3Max, a4, b4);
a4 = _mm_shufflelo_epi16(a4, _MM_SHUFFLE(0, 1, 2, 3));
b4 = _mm_shufflehi_epi16(b4, _MM_SHUFFLE(0, 1, 2, 3));
__m128i a4Min = _mm_min_epi16(a4, b4);
__m128i b4Max = _mm_max_epi16(a4, b4);
__m128i a5, b5;
unpack_epi16_ab(a4Min, b4Max, a5, b5);
unpack_epi64_ab(a5, b5, a5, b5);
a5 = _mm_shuffle_epi32(a5, _MM_SHUFFLE(2, 3, 0, 1));
__m128i a5Min = _mm_min_epi16(a5, b5);
__m128i b5Max = _mm_max_epi16(a5, b5);
__m128i a6, b6;
unpack_epi16_ab(a5Min, b5Max, a6, b6);
unpack_epi16_ab(a6, b6, a6, b6);
__m128i a6Min = _mm_min_epi16(a6, b6);
__m128i b6Max = _mm_max_epi16(a6, b6);
__m128i a7, b7;
b7 = _mm_shuffle_epi32(b6Max, _MM_SHUFFLE(1, 0, 3, 2));
b7 = _mm_shufflehi_epi16(b7, _MM_SHUFFLE(0, 1, 2, 3));
a7 = _mm_shufflelo_epi16(a6Min, _MM_SHUFFLE(0, 1, 2, 3));
__m128i a7Min = _mm_min_epi16(a7, b7);
__m128i b7Max = _mm_max_epi16(a7, b7);
__m128i a8, b8;
unpack_epi16_ab(a7Min, b7Max, b8, a8);
a8 = _mm_shuffle_epi32(a8, _MM_SHUFFLE(0, 1, 2, 3));
__m128i a8Min = _mm_min_epi16(a8, b8);
__m128i b8Max = _mm_max_epi16(a8, b8);
__m128i a9, b9;
unpack_epi16_ab(a8Min, b8Max, a9, b9);
__m128i a9Min = _mm_min_epi16(a9, b9);
__m128i b9Max = _mm_max_epi16(a9, b9);
__m128i a10, b10;
unpack_epi16_ab(a9Min, b9Max, a10, b10);
__m128i a10Min = _mm_min_epi16(a10, b10);
__m128i b10Max = _mm_max_epi16(a10, b10);
unpack_epi16_ab(a10Min, b10Max, aOut, bOut);
}
This is verified using the standard 0/1 technique described in Knuth.
The compiled code uses only five registers, so avoids register spills.
I am still not sure if this was worth the effort. But it's done.
Sunday, March 19, 2017
shaving a few shuffles from the float sorting network in SSE2
In an earlier post I commented that I felt that it should be possible to get rid of a couple shuffles in the presented network.
Below is a network implementation of the same bitonic network above. But this new implementation seems to have a handful fewer shuffles.
void BitonicSortingNetwork_FewerShuffles_SSE2(const __m128 &a, const __m128 &b, __m128 &aOut, __m128 &bOut)
{
__m128 a1 = a;
__m128 b1 = b;
__m128 a1Min = _mm_min_ps(a1, b1); // 1357
__m128 b1Max = _mm_max_ps(a1, b1); // 2468
__m128 a2 = a1Min; // 1357
__m128 b2 = _mm_shuffle_ps(b1Max, b1Max, _MM_SHUFFLE(2, 3, 0, 1)); // 4286
__m128 a2Min = _mm_min_ps(a2, b2); // 1256
__m128 b2Max = _mm_max_ps(a2, b2); // 4387
__m128 a3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(3, 1, 2, 0)); // 1537
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(2, 0, 3,1)); // 2648
__m128 a3Min = _mm_min_ps(a3, b3); // 1537
__m128 b3Max = _mm_max_ps(a3, b3); // 2648
__m128 a4 = a3Min; // 1537
__m128 b4 = _mm_shuffle_ps(b3Max, b3Max, _MM_SHUFFLE(0, 1, 2, 3)); // 8462
__m128 a4Min = _mm_min_ps(a4, b4); // 1432
__m128 b4Max = _mm_max_ps(a4, b4); // 8567
__m128 a5 = _mm_unpacklo_ps(a4Min, b4Max); // 1845
__m128 b5 = _mm_unpackhi_ps(a4Min, b4Max); // 3627
__m128 a5Min = _mm_min_ps(a5, b5); // 1625
__m128 b5Max = _mm_max_ps(a5, b5); // 3847
__m128 a6 = _mm_unpacklo_ps(a5Min, b5Max); // 1368
__m128 b6 = _mm_unpackhi_ps(a5Min, b5Max); // 2457
__m128 a6Min = _mm_min_ps(a6, b6); // 1357
__m128 b6Max = _mm_max_ps(a6, b6); // 2468
aOut = _mm_unpacklo_ps(a6Min, b6Max); // 1234
bOut = _mm_unpackhi_ps(a6Min, b6Max); // 5678
}
In particular it arranges to have only the unpacks before returning the results. By re-using an intermediate, two shuffles are eliminated between min/max. While this only shaves a cycle or two from the runtime, it is more aesthetically pleasing.
The compiled assembly code looks good.
__m128 a1 = a;
00007FF751C51070 movups xmm3,xmmword ptr [rcx]
__m128 b1 = b;
__m128 a1Min = _mm_min_ps(a1, b1); // 1357
__m128 b1Max = _mm_max_ps(a1, b1); // 2468
00007FF751C51073 movaps xmm1,xmm3
00007FF751C51076 maxps xmm1,xmmword ptr [rdx]
00007FF751C51079 minps xmm3,xmmword ptr [rdx]
__m128 a2 = a1Min; // 1357
__m128 b2 = _mm_shuffle_ps(b1Max, b1Max, _MM_SHUFFLE(2, 3, 0, 1)); // 4286
00007FF751C5107C shufps xmm1,xmm1,0B1h
__m128 a2Min = _mm_min_ps(a2, b2); // 1256
00007FF751C51080 movaps xmm2,xmm3
00007FF751C51083 minps xmm2,xmm1
__m128 b2Max = _mm_max_ps(a2, b2); // 4387
00007FF751C51086 maxps xmm3,xmm1
__m128 a3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(3, 1, 2, 0)); // 1537
00007FF751C51089 movaps xmm4,xmm2
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(2, 0, 3, 1)); // 2648
00007FF751C5108C shufps xmm2,xmm3,8Dh
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(2, 0, 3, 1)); // 2648
00007FF751C51090 shufps xmm4,xmm3,0D8h
__m128 a3Min = _mm_min_ps(a3, b3); // 1537
__m128 b3Max = _mm_max_ps(a3, b3); // 2648
00007FF751C51094 movaps xmm0,xmm4
00007FF751C51097 maxps xmm0,xmm2
00007FF751C5109A minps xmm4,xmm2
__m128 a4 = a3Min; // 1537
__m128 b4 = _mm_shuffle_ps(b3Max, b3Max, _MM_SHUFFLE(0, 1, 2, 3)); // 8462
00007FF751C5109D shufps xmm0,xmm0,1Bh
__m128 a4Min = _mm_min_ps(a4, b4); // 1432
00007FF751C510A1 movaps xmm1,xmm4
00007FF751C510A4 minps xmm1,xmm0
__m128 b4Max = _mm_max_ps(a4, b4); // 8567
00007FF751C510A7 maxps xmm4,xmm0
__m128 a5 = _mm_unpacklo_ps(a4Min, b4Max); // 1845
00007FF751C510AA movaps xmm0,xmm1
__m128 b5 = _mm_unpackhi_ps(a4Min, b4Max); // 3627
00007FF751C510AD unpckhps xmm1,xmm4
00007FF751C510B0 unpcklps xmm0,xmm4
__m128 a5Min = _mm_min_ps(a5, b5); // 1625
00007FF751C510B3 movaps xmm2,xmm0
__m128 b5Max = _mm_max_ps(a5, b5); // 3847
00007FF751C510B6 maxps xmm0,xmm1
00007FF751C510B9 minps xmm2,xmm1
__m128 a6 = _mm_unpacklo_ps(a5Min, b5Max); // 1368
00007FF751C510BC movaps xmm3,xmm2
__m128 b6 = _mm_unpackhi_ps(a5Min, b5Max); // 2457
00007FF751C510BF unpckhps xmm2,xmm0
00007FF751C510C2 unpcklps xmm3,xmm0
__m128 a6Min = _mm_min_ps(a6, b6); // 1357
00007FF751C510C5 movaps xmm1,xmm3
00007FF751C510C8 minps xmm1,xmm2
__m128 b6Max = _mm_max_ps(a6, b6); // 2468
00007FF751C510CB maxps xmm3,xmm2
aOut = _mm_unpacklo_ps(a6Min, b6Max); // 1234
00007FF751C510CE movaps xmm0,xmm1
00007FF751C510D1 unpcklps xmm0,xmm3
bOut = _mm_unpackhi_ps(a6Min, b6Max); // 5678
00007FF751C510D4 unpckhps xmm1,xmm3
00007FF751C510D7 movups xmmword ptr [r8],xmm0
00007FF751C510DB movups xmmword ptr [r9],xmm1
}
00007FF751C510DF ret
Update: I see that the final three stages are the same as what is described in a paper that I read but hadn't remembered. I find a different paper in a quick search. I think that the first three stages are better than what are described in that paper. But incorporated into a much larger mergesort gives a different set of problems. Also, the papers don't give source code...
Below is a network implementation of the same bitonic network above. But this new implementation seems to have a handful fewer shuffles.
void BitonicSortingNetwork_FewerShuffles_SSE2(const __m128 &a, const __m128 &b, __m128 &aOut, __m128 &bOut)
{
__m128 a1 = a;
__m128 b1 = b;
__m128 a1Min = _mm_min_ps(a1, b1); // 1357
__m128 b1Max = _mm_max_ps(a1, b1); // 2468
__m128 a2 = a1Min; // 1357
__m128 b2 = _mm_shuffle_ps(b1Max, b1Max, _MM_SHUFFLE(2, 3, 0, 1)); // 4286
__m128 a2Min = _mm_min_ps(a2, b2); // 1256
__m128 b2Max = _mm_max_ps(a2, b2); // 4387
__m128 a3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(3, 1, 2, 0)); // 1537
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(2, 0, 3,1)); // 2648
__m128 a3Min = _mm_min_ps(a3, b3); // 1537
__m128 b3Max = _mm_max_ps(a3, b3); // 2648
__m128 a4 = a3Min; // 1537
__m128 b4 = _mm_shuffle_ps(b3Max, b3Max, _MM_SHUFFLE(0, 1, 2, 3)); // 8462
__m128 a4Min = _mm_min_ps(a4, b4); // 1432
__m128 b4Max = _mm_max_ps(a4, b4); // 8567
__m128 a5 = _mm_unpacklo_ps(a4Min, b4Max); // 1845
__m128 b5 = _mm_unpackhi_ps(a4Min, b4Max); // 3627
__m128 a5Min = _mm_min_ps(a5, b5); // 1625
__m128 b5Max = _mm_max_ps(a5, b5); // 3847
__m128 a6 = _mm_unpacklo_ps(a5Min, b5Max); // 1368
__m128 b6 = _mm_unpackhi_ps(a5Min, b5Max); // 2457
__m128 a6Min = _mm_min_ps(a6, b6); // 1357
__m128 b6Max = _mm_max_ps(a6, b6); // 2468
aOut = _mm_unpacklo_ps(a6Min, b6Max); // 1234
bOut = _mm_unpackhi_ps(a6Min, b6Max); // 5678
}
In particular it arranges to have only the unpacks before returning the results. By re-using an intermediate, two shuffles are eliminated between min/max. While this only shaves a cycle or two from the runtime, it is more aesthetically pleasing.
The compiled assembly code looks good.
__m128 a1 = a;
00007FF751C51070 movups xmm3,xmmword ptr [rcx]
__m128 b1 = b;
__m128 a1Min = _mm_min_ps(a1, b1); // 1357
__m128 b1Max = _mm_max_ps(a1, b1); // 2468
00007FF751C51073 movaps xmm1,xmm3
00007FF751C51076 maxps xmm1,xmmword ptr [rdx]
00007FF751C51079 minps xmm3,xmmword ptr [rdx]
__m128 a2 = a1Min; // 1357
__m128 b2 = _mm_shuffle_ps(b1Max, b1Max, _MM_SHUFFLE(2, 3, 0, 1)); // 4286
00007FF751C5107C shufps xmm1,xmm1,0B1h
__m128 a2Min = _mm_min_ps(a2, b2); // 1256
00007FF751C51080 movaps xmm2,xmm3
00007FF751C51083 minps xmm2,xmm1
__m128 b2Max = _mm_max_ps(a2, b2); // 4387
00007FF751C51086 maxps xmm3,xmm1
__m128 a3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(3, 1, 2, 0)); // 1537
00007FF751C51089 movaps xmm4,xmm2
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(2, 0, 3, 1)); // 2648
00007FF751C5108C shufps xmm2,xmm3,8Dh
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(2, 0, 3, 1)); // 2648
00007FF751C51090 shufps xmm4,xmm3,0D8h
__m128 a3Min = _mm_min_ps(a3, b3); // 1537
__m128 b3Max = _mm_max_ps(a3, b3); // 2648
00007FF751C51094 movaps xmm0,xmm4
00007FF751C51097 maxps xmm0,xmm2
00007FF751C5109A minps xmm4,xmm2
__m128 a4 = a3Min; // 1537
__m128 b4 = _mm_shuffle_ps(b3Max, b3Max, _MM_SHUFFLE(0, 1, 2, 3)); // 8462
00007FF751C5109D shufps xmm0,xmm0,1Bh
__m128 a4Min = _mm_min_ps(a4, b4); // 1432
00007FF751C510A1 movaps xmm1,xmm4
00007FF751C510A4 minps xmm1,xmm0
__m128 b4Max = _mm_max_ps(a4, b4); // 8567
00007FF751C510A7 maxps xmm4,xmm0
__m128 a5 = _mm_unpacklo_ps(a4Min, b4Max); // 1845
00007FF751C510AA movaps xmm0,xmm1
__m128 b5 = _mm_unpackhi_ps(a4Min, b4Max); // 3627
00007FF751C510AD unpckhps xmm1,xmm4
00007FF751C510B0 unpcklps xmm0,xmm4
__m128 a5Min = _mm_min_ps(a5, b5); // 1625
00007FF751C510B3 movaps xmm2,xmm0
__m128 b5Max = _mm_max_ps(a5, b5); // 3847
00007FF751C510B6 maxps xmm0,xmm1
00007FF751C510B9 minps xmm2,xmm1
__m128 a6 = _mm_unpacklo_ps(a5Min, b5Max); // 1368
00007FF751C510BC movaps xmm3,xmm2
__m128 b6 = _mm_unpackhi_ps(a5Min, b5Max); // 2457
00007FF751C510BF unpckhps xmm2,xmm0
00007FF751C510C2 unpcklps xmm3,xmm0
__m128 a6Min = _mm_min_ps(a6, b6); // 1357
00007FF751C510C5 movaps xmm1,xmm3
00007FF751C510C8 minps xmm1,xmm2
__m128 b6Max = _mm_max_ps(a6, b6); // 2468
00007FF751C510CB maxps xmm3,xmm2
aOut = _mm_unpacklo_ps(a6Min, b6Max); // 1234
00007FF751C510CE movaps xmm0,xmm1
00007FF751C510D1 unpcklps xmm0,xmm3
bOut = _mm_unpackhi_ps(a6Min, b6Max); // 5678
00007FF751C510D4 unpckhps xmm1,xmm3
00007FF751C510D7 movups xmmword ptr [r8],xmm0
00007FF751C510DB movups xmmword ptr [r9],xmm1
}
00007FF751C510DF ret
Wednesday, March 15, 2017
A bitonic sorting network for 8 floats in SSE2
The previous sorting network did perform the sort in the ideal number of min/max operations. But there was a lot of shuffling. But this did strict ordering of the input wires. If instead we allow the input wires to be reordered to make the shuffles more efficient, we can do a little better.
The below is a bitonic sorting network performed in this manner. We shuffle to get the correct comparisons without regard to where in the registers the conceptual input wires are located. This verifies as correct with the standard 0/1 method.
The below network is done with floats, since all of the necessary intrisics are available in SSE2.
void BitonicSortingNetwork_SSE2(const __m128 &a, const __m128 &b, __m128 &aOut, __m128 &bOut)
{
__m128 a1 = _mm_shuffle_ps(a, b, _MM_SHUFFLE(2, 0, 2, 0)); // 1357
__m128 b1 = _mm_shuffle_ps(a, b, _MM_SHUFFLE(3, 1, 3, 1)); // 2468
__m128 a1Min = _mm_min_ps(a1, b1);
__m128 b1Max = _mm_max_ps(a1, b1);
__m128 a2 = _mm_shuffle_ps(a1Min, b1Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1526
__m128 b2 = _mm_shuffle_ps(b1Max, a1Min, _MM_SHUFFLE(3, 1, 3, 1)); // 4837
__m128 a2Min = _mm_min_ps(a2, b2);
__m128 b2Max = _mm_max_ps(a2, b2);
__m128 a3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(3, 2, 1, 0)); // 1537
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(1, 0, 3, 2)); // 2648
__m128 a3Min = _mm_min_ps(a3, b3);
__m128 b3Max = _mm_max_ps(a3, b3);
__m128 a4 = _mm_shuffle_ps(a3Min, b3Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1324
__m128 b4 = _mm_shuffle_ps(b3Max, a3Min, _MM_SHUFFLE(1, 3, 1, 3)); // 8675
__m128 a4Min = _mm_min_ps(a4, b4);
__m128 b4Max = _mm_max_ps(a4, b4);
__m128 a5 = _mm_shuffle_ps(a4Min, b4Max, _MM_SHUFFLE(1, 3, 2, 0)); // 1256
__m128 b5 = _mm_shuffle_ps(a4Min, b4Max, _MM_SHUFFLE(0, 2, 3, 1)); // 3478
__m128 a5Min = _mm_min_ps(a5, b5);
__m128 b5Max = _mm_max_ps(a5, b5);
__m128 a6 = _mm_shuffle_ps(a5Min, b5Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1537
__m128 b6 = _mm_shuffle_ps(a5Min, b5Max, _MM_SHUFFLE(3, 1, 3, 1)); // 2648
__m128 a6Min = _mm_min_ps(a6, b6);
__m128 b6Max = _mm_max_ps(a6, b6);
__m128 a7 = _mm_shuffle_ps(a6Min, b6Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1324
__m128 b7 = _mm_shuffle_ps(a6Min, b6Max, _MM_SHUFFLE(3, 1, 3, 1)); // 5768
aOut = _mm_shuffle_ps(a7, a7, _MM_SHUFFLE(3, 1, 2, 0)); // 1234
bOut = _mm_shuffle_ps(b7, b7, _MM_SHUFFLE(3, 1, 2, 0)); // 5678
}
The comments reflect where in the registers the conceptual input wires are located.
Since this is in fact a sorting network, the order of the inputs isn't important, and the first pair of shuffles is unnecessary. I've tried to figure out a way of getting rid of the final double shuffle, but I haven't been able to get rid of it.
Here's the assembly of the compilation without the initial shuffles:
__m128 a1 = a;
000000013F1B1140 movaps xmm2,xmmword ptr [rcx]
__m128 b1 = b;
__m128 a1Min = _mm_min_ps(a1, b1);
000000013F1B1143 movaps xmm1,xmm2
__m128 b1Max = _mm_max_ps(a1, b1);
000000013F1B1146 maxps xmm2,xmmword ptr [rdx]
000000013F1B1149 minps xmm1,xmmword ptr [rdx]
__m128 a2 = _mm_shuffle_ps(a1Min, b1Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1526
000000013F1B114C movaps xmm0,xmm1
000000013F1B114F shufps xmm0,xmm2,88h
__m128 b2 = _mm_shuffle_ps(b1Max, a1Min, _MM_SHUFFLE(3, 1, 3, 1)); // 4837
000000013F1B1153 shufps xmm2,xmm1,0DDh
__m128 a2Min = _mm_min_ps(a2, b2);
000000013F1B1157 movaps xmm1,xmm0
__m128 b2Max = _mm_max_ps(a2, b2);
000000013F1B115A maxps xmm0,xmm2
000000013F1B115D minps xmm1,xmm2
__m128 a3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(3, 2, 1, 0)); // 1537
000000013F1B1160 movaps xmm2,xmm1
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(1, 0, 3, 2)); // 2648
000000013F1B1163 shufps xmm1,xmm0,4Eh
000000013F1B1167 shufps xmm2,xmm0,0E4h
__m128 a3Min = _mm_min_ps(a3, b3);
000000013F1B116B movaps xmm0,xmm2
__m128 b3Max = _mm_max_ps(a3, b3);
000000013F1B116E maxps xmm2,xmm1
000000013F1B1171 minps xmm0,xmm1
__m128 a4 = _mm_shuffle_ps(a3Min, b3Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1324
000000013F1B1174 movaps xmm1,xmm0
000000013F1B1177 shufps xmm1,xmm2,88h
__m128 b4 = _mm_shuffle_ps(b3Max, a3Min, _MM_SHUFFLE(1, 3, 1, 3)); // 8675
000000013F1B117B shufps xmm2,xmm0,77h
__m128 a4Min = _mm_min_ps(a4, b4);
000000013F1B117F movaps xmm0,xmm1
000000013F1B1182 minps xmm0,xmm2
__m128 b4Max = _mm_max_ps(a4, b4);
000000013F1B1185 maxps xmm1,xmm2
__m128 a5 = _mm_shuffle_ps(a4Min, b4Max, _MM_SHUFFLE(1, 3, 2, 0)); // 1256
000000013F1B1188 movaps xmm2,xmm0
000000013F1B118B shufps xmm2,xmm1,78h
__m128 b5 = _mm_shuffle_ps(a4Min, b4Max, _MM_SHUFFLE(0, 2, 3, 1)); // 3478
000000013F1B118F shufps xmm0,xmm1,2Dh
__m128 a5Min = _mm_min_ps(a5, b5);
000000013F1B1193 movaps xmm1,xmm2
__m128 b5Max = _mm_max_ps(a5, b5);
000000013F1B1196 maxps xmm2,xmm0
000000013F1B1199 minps xmm1,xmm0
__m128 a6 = _mm_shuffle_ps(a5Min, b5Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1537
000000013F1B119C movaps xmm3,xmm1
__m128 b6 = _mm_shuffle_ps(a5Min, b5Max, _MM_SHUFFLE(3, 1, 3, 1)); // 2648
000000013F1B119F shufps xmm1,xmm2,0DDh
000000013F1B11A3 shufps xmm3,xmm2,88h
__m128 a6Min = _mm_min_ps(a6, b6);
000000013F1B11A7 movaps xmm2,xmm3
000000013F1B11AA minps xmm2,xmm1
__m128 b6Max = _mm_max_ps(a6, b6);
000000013F1B11AD maxps xmm3,xmm1
__m128 a7 = _mm_shuffle_ps(a6Min, b6Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1324
000000013F1B11B0 movaps xmm0,xmm2
000000013F1B11B3 shufps xmm0,xmm3,88h
__m128 b7 = _mm_shuffle_ps(a6Min, b6Max, _MM_SHUFFLE(3, 1, 3, 1)); // 5768
000000013F1B11B7 shufps xmm2,xmm3,0DDh
aOut = _mm_shuffle_ps(a7, a7, _MM_SHUFFLE(3, 1, 2, 0)); // 1234
000000013F1B11BB shufps xmm0,xmm0,0D8h
bOut = _mm_shuffle_ps(b7, b7, _MM_SHUFFLE(3, 1, 2, 0)); // 5678
000000013F1B11BF shufps xmm2,xmm2,0D8h
000000013F1B11C3 movaps xmmword ptr [r8],xmm0
000000013F1B11C7 movaps xmmword ptr [r9],xmm2
}
000000013F1B11CB ret
This fits snugly into 4 registers. It looks pretty good.
The below is a bitonic sorting network performed in this manner. We shuffle to get the correct comparisons without regard to where in the registers the conceptual input wires are located. This verifies as correct with the standard 0/1 method.
The below network is done with floats, since all of the necessary intrisics are available in SSE2.
void BitonicSortingNetwork_SSE2(const __m128 &a, const __m128 &b, __m128 &aOut, __m128 &bOut)
{
__m128 a1 = _mm_shuffle_ps(a, b, _MM_SHUFFLE(2, 0, 2, 0)); // 1357
__m128 b1 = _mm_shuffle_ps(a, b, _MM_SHUFFLE(3, 1, 3, 1)); // 2468
__m128 a1Min = _mm_min_ps(a1, b1);
__m128 b1Max = _mm_max_ps(a1, b1);
__m128 a2 = _mm_shuffle_ps(a1Min, b1Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1526
__m128 b2 = _mm_shuffle_ps(b1Max, a1Min, _MM_SHUFFLE(3, 1, 3, 1)); // 4837
__m128 a2Min = _mm_min_ps(a2, b2);
__m128 b2Max = _mm_max_ps(a2, b2);
__m128 a3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(3, 2, 1, 0)); // 1537
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(1, 0, 3, 2)); // 2648
__m128 a3Min = _mm_min_ps(a3, b3);
__m128 b3Max = _mm_max_ps(a3, b3);
__m128 a4 = _mm_shuffle_ps(a3Min, b3Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1324
__m128 b4 = _mm_shuffle_ps(b3Max, a3Min, _MM_SHUFFLE(1, 3, 1, 3)); // 8675
__m128 a4Min = _mm_min_ps(a4, b4);
__m128 b4Max = _mm_max_ps(a4, b4);
__m128 a5 = _mm_shuffle_ps(a4Min, b4Max, _MM_SHUFFLE(1, 3, 2, 0)); // 1256
__m128 b5 = _mm_shuffle_ps(a4Min, b4Max, _MM_SHUFFLE(0, 2, 3, 1)); // 3478
__m128 a5Min = _mm_min_ps(a5, b5);
__m128 b5Max = _mm_max_ps(a5, b5);
__m128 a6 = _mm_shuffle_ps(a5Min, b5Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1537
__m128 b6 = _mm_shuffle_ps(a5Min, b5Max, _MM_SHUFFLE(3, 1, 3, 1)); // 2648
__m128 a6Min = _mm_min_ps(a6, b6);
__m128 b6Max = _mm_max_ps(a6, b6);
__m128 a7 = _mm_shuffle_ps(a6Min, b6Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1324
__m128 b7 = _mm_shuffle_ps(a6Min, b6Max, _MM_SHUFFLE(3, 1, 3, 1)); // 5768
aOut = _mm_shuffle_ps(a7, a7, _MM_SHUFFLE(3, 1, 2, 0)); // 1234
bOut = _mm_shuffle_ps(b7, b7, _MM_SHUFFLE(3, 1, 2, 0)); // 5678
}
The comments reflect where in the registers the conceptual input wires are located.
Since this is in fact a sorting network, the order of the inputs isn't important, and the first pair of shuffles is unnecessary. I've tried to figure out a way of getting rid of the final double shuffle, but I haven't been able to get rid of it.
Here's the assembly of the compilation without the initial shuffles:
__m128 a1 = a;
000000013F1B1140 movaps xmm2,xmmword ptr [rcx]
__m128 b1 = b;
__m128 a1Min = _mm_min_ps(a1, b1);
000000013F1B1143 movaps xmm1,xmm2
__m128 b1Max = _mm_max_ps(a1, b1);
000000013F1B1146 maxps xmm2,xmmword ptr [rdx]
000000013F1B1149 minps xmm1,xmmword ptr [rdx]
__m128 a2 = _mm_shuffle_ps(a1Min, b1Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1526
000000013F1B114C movaps xmm0,xmm1
000000013F1B114F shufps xmm0,xmm2,88h
__m128 b2 = _mm_shuffle_ps(b1Max, a1Min, _MM_SHUFFLE(3, 1, 3, 1)); // 4837
000000013F1B1153 shufps xmm2,xmm1,0DDh
__m128 a2Min = _mm_min_ps(a2, b2);
000000013F1B1157 movaps xmm1,xmm0
__m128 b2Max = _mm_max_ps(a2, b2);
000000013F1B115A maxps xmm0,xmm2
000000013F1B115D minps xmm1,xmm2
__m128 a3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(3, 2, 1, 0)); // 1537
000000013F1B1160 movaps xmm2,xmm1
__m128 b3 = _mm_shuffle_ps(a2Min, b2Max, _MM_SHUFFLE(1, 0, 3, 2)); // 2648
000000013F1B1163 shufps xmm1,xmm0,4Eh
000000013F1B1167 shufps xmm2,xmm0,0E4h
__m128 a3Min = _mm_min_ps(a3, b3);
000000013F1B116B movaps xmm0,xmm2
__m128 b3Max = _mm_max_ps(a3, b3);
000000013F1B116E maxps xmm2,xmm1
000000013F1B1171 minps xmm0,xmm1
__m128 a4 = _mm_shuffle_ps(a3Min, b3Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1324
000000013F1B1174 movaps xmm1,xmm0
000000013F1B1177 shufps xmm1,xmm2,88h
__m128 b4 = _mm_shuffle_ps(b3Max, a3Min, _MM_SHUFFLE(1, 3, 1, 3)); // 8675
000000013F1B117B shufps xmm2,xmm0,77h
__m128 a4Min = _mm_min_ps(a4, b4);
000000013F1B117F movaps xmm0,xmm1
000000013F1B1182 minps xmm0,xmm2
__m128 b4Max = _mm_max_ps(a4, b4);
000000013F1B1185 maxps xmm1,xmm2
__m128 a5 = _mm_shuffle_ps(a4Min, b4Max, _MM_SHUFFLE(1, 3, 2, 0)); // 1256
000000013F1B1188 movaps xmm2,xmm0
000000013F1B118B shufps xmm2,xmm1,78h
__m128 b5 = _mm_shuffle_ps(a4Min, b4Max, _MM_SHUFFLE(0, 2, 3, 1)); // 3478
000000013F1B118F shufps xmm0,xmm1,2Dh
__m128 a5Min = _mm_min_ps(a5, b5);
000000013F1B1193 movaps xmm1,xmm2
__m128 b5Max = _mm_max_ps(a5, b5);
000000013F1B1196 maxps xmm2,xmm0
000000013F1B1199 minps xmm1,xmm0
__m128 a6 = _mm_shuffle_ps(a5Min, b5Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1537
000000013F1B119C movaps xmm3,xmm1
__m128 b6 = _mm_shuffle_ps(a5Min, b5Max, _MM_SHUFFLE(3, 1, 3, 1)); // 2648
000000013F1B119F shufps xmm1,xmm2,0DDh
000000013F1B11A3 shufps xmm3,xmm2,88h
__m128 a6Min = _mm_min_ps(a6, b6);
000000013F1B11A7 movaps xmm2,xmm3
000000013F1B11AA minps xmm2,xmm1
__m128 b6Max = _mm_max_ps(a6, b6);
000000013F1B11AD maxps xmm3,xmm1
__m128 a7 = _mm_shuffle_ps(a6Min, b6Max, _MM_SHUFFLE(2, 0, 2, 0)); // 1324
000000013F1B11B0 movaps xmm0,xmm2
000000013F1B11B3 shufps xmm0,xmm3,88h
__m128 b7 = _mm_shuffle_ps(a6Min, b6Max, _MM_SHUFFLE(3, 1, 3, 1)); // 5768
000000013F1B11B7 shufps xmm2,xmm3,0DDh
aOut = _mm_shuffle_ps(a7, a7, _MM_SHUFFLE(3, 1, 2, 0)); // 1234
000000013F1B11BB shufps xmm0,xmm0,0D8h
bOut = _mm_shuffle_ps(b7, b7, _MM_SHUFFLE(3, 1, 2, 0)); // 5678
000000013F1B11BF shufps xmm2,xmm2,0D8h
000000013F1B11C3 movaps xmmword ptr [r8],xmm0
000000013F1B11C7 movaps xmmword ptr [r9],xmm2
}
000000013F1B11CB ret
This fits snugly into 4 registers. It looks pretty good.
Wednesday, March 1, 2017
Perfect shuffles of 16bit numbers in two registers
In a single 128bit register, one can store eight 16bit numbers. For a pair of these 128bit registers, the instruction pair _mm_unpacklo_epi16 and _mm_unpackhi_epi16, called with the same two registers as arguments, will conceptually perform a perfect shuffle of the 16-bit numbers in the registers.
void printm_epu16(const __m128i &mx)
{
for (int iA = 0; iA < ARRAYSIZE(mx.m128i_u16); iA++)
{
printf("%x", mx.m128i_u16[iA]);
}
}
int main()
{
__m128i a = _mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7);
__m128i b = _mm_setr_epi16(8, 9, 10, 11, 12, 13, 14, 15);
printm_epu16(a);
printm_epu16(b);
printf("\n");
for (int i = 0; i < 4; i++)
{
__m128i aa = _mm_unpacklo_epi16(a, b);
__m128i bb = _mm_unpackhi_epi16(a, b);
a = aa;
b = bb;
printm_epu16(a);
printm_epu16(b);
printf("\n");
}
}
will yield the output:
0123456789abcdef
08192a3b4c5d6e7f
048c159d26ae37bf
02468ace13579bdf
0123456789abcdef
We see that this cycle repeats after four iterations.
void printm_epu16(const __m128i &mx)
{
for (int iA = 0; iA < ARRAYSIZE(mx.m128i_u16); iA++)
{
printf("%x", mx.m128i_u16[iA]);
}
}
int main()
{
__m128i a = _mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7);
__m128i b = _mm_setr_epi16(8, 9, 10, 11, 12, 13, 14, 15);
printm_epu16(a);
printm_epu16(b);
printf("\n");
for (int i = 0; i < 4; i++)
{
__m128i aa = _mm_unpacklo_epi16(a, b);
__m128i bb = _mm_unpackhi_epi16(a, b);
a = aa;
b = bb;
printm_epu16(a);
printm_epu16(b);
printf("\n");
}
}
will yield the output:
0123456789abcdef
08192a3b4c5d6e7f
048c159d26ae37bf
02468ace13579bdf
0123456789abcdef
We see that this cycle repeats after four iterations.
Monday, January 16, 2017
transposing a non-square float matrix with SSE
The transpose of a matrix can be useful. During a matrix multiply a transpose of the right hand matrix will allow the result to be expressed as a series of dot products, with the inputs of the dot product sequential in memory.
But often the matrix will be non-square, and so cannot be transposed in place. But let's look at transposing with a source and destination.
We will transpose in 4x4 blocks. Momentarily ignoring what is evenly divisible by 4:
for (size_t y = nHeight&~3; y >0;)
{
y -= 4;
for (size_t x = nWidth&~3; x>0;)
{
x -= 4;
__m128 row0;
__m128 row1;
__m128 row2;
__m128 row3;
row0 = _mm_loadu_ps(&pMatrixIn[(y + 0)*nWidth + x]);
row1 = _mm_loadu_ps(&pMatrixIn[(y + 1)*nWidth + x]);
row2 = _mm_loadu_ps(&pMatrixIn[(y + 2)*nWidth + x]);
row3 = _mm_loadu_ps(&pMatrixIn[(y + 3)*nWidth + x]);
_MM_TRANSPOSE4_PS(row0, row1, row2, row3);
_mm_storeu_ps(&pMatrixOut[(x + 0)*nHeight + y], row0);
_mm_storeu_ps(&pMatrixOut[(x + 1)*nHeight + y], row1);
_mm_storeu_ps(&pMatrixOut[(x + 2)*nHeight + y], row2);
_mm_storeu_ps(&pMatrixOut[(x + 3)*nHeight + y], row3);
}
}
This seems pretty tight...but that's a lot of multiplies. We can factor out most of the multiplies so that it's just pointer arithmetic and adds. On x64 the central loop can be accomplished without multiplies and without register spills.
Most of the matrix will be transposed by this set of loops. We will have to handle the rest in two other rectangular strips (making sure not to double-count the corner rectangle).
void TransposeNonSquareMatrix_SSE(const float *const pMatrixIn, float *const pMatrixOut, const size_t nWidth, const size_t nHeight)
{
for (size_t y = nHeight; y & 3;)
{
y--;
for (size_t x = nWidth; x-- > 0;)
{
pMatrixOut[x*nHeight + y] = pMatrixIn[y*nWidth + x];
}
}
if (nWidth & 3)
{
for (size_t y = nHeight&~3; y-- > 0;)
{
for (size_t x = nWidth; x & 3;)
{
x--;
pMatrixOut[x*nHeight + y] = pMatrixIn[y*nWidth + x];
}
}
}
const size_t nHeightMax = nHeight&~3;
const size_t nWidthMax = nWidth&~3;
for (size_t y = nHeightMax; y >0;)
{
y -= 4;
const float *const pIn0 = pMatrixIn + y*nWidth;
const float *const pIn1 = pIn0 + nWidth;
const float *const pIn2 = pIn1 + nWidth;
const float * const pIn3 = pIn2 + nWidth;
float * pOut = &pMatrixOut[nWidthMax*nHeight + y];
for (size_t x = nWidthMax; x>0;)
{
x -= 4;
__m128 row0;
__m128 row1;
__m128 row2;
__m128 row3;
row0 = _mm_loadu_ps(&pIn0[x]);
row1 = _mm_loadu_ps(&pIn1[x]);
row2 = _mm_loadu_ps(&pIn2[x]);
row3 = _mm_loadu_ps(&pIn3[x]);
_MM_TRANSPOSE4_PS(row0, row1, row2, row3);
_mm_storeu_ps((pOut -= nHeight), row3);
_mm_storeu_ps((pOut -= nHeight), row2);
_mm_storeu_ps((pOut -= nHeight), row1);
_mm_storeu_ps((pOut -= nHeight), row0);
}
}
}
The assembly compilation of the central loop is pretty tight:
{
x -= 4;
__m128 row0;
__m128 row1;
__m128 row2;
__m128 row3;
row0 = _mm_loadu_ps(&pIn0[x]);
row1 = _mm_loadu_ps(&pIn1[x]);
000000013F401210 movups xmm0,xmmword ptr [r10-10h]
row2 = _mm_loadu_ps(&pIn2[x]);
row3 = _mm_loadu_ps(&pIn3[x]);
_MM_TRANSPOSE4_PS(row0, row1, row2, row3);
_mm_storeu_ps((pOut -= nHeight), row3);
000000013F401215 sub rdx,rbp
000000013F401218 lea r10,[r10-10h]
000000013F40121C movups xmm4,xmmword ptr [r10+r15]
000000013F401221 movups xmm3,xmmword ptr [r10+r11]
000000013F401226 movups xmm1,xmmword ptr [r10+rcx]
000000013F40122B movaps xmm5,xmm4
000000013F40122E movaps xmm2,xmm3
000000013F401231 shufps xmm4,xmm0,0EEh
000000013F401235 shufps xmm5,xmm0,44h
row2 = _mm_loadu_ps(&pIn2[x]);
row3 = _mm_loadu_ps(&pIn3[x]);
_MM_TRANSPOSE4_PS(row0, row1, row2, row3);
_mm_storeu_ps((pOut -= nHeight), row3);
000000013F401239 movaps xmm0,xmm4
000000013F40123C shufps xmm3,xmm1,0EEh
000000013F401240 shufps xmm0,xmm3,0DDh
000000013F401244 movups xmmword ptr [rdx],xmm0
_mm_storeu_ps((pOut -= nHeight), row2);
000000013F401247 sub rdx,rbp
000000013F40124A shufps xmm2,xmm1,44h
000000013F40124E movaps xmm0,xmm5
000000013F401251 shufps xmm4,xmm3,88h
000000013F401255 shufps xmm0,xmm2,0DDh
000000013F401259 movups xmmword ptr [rdx],xmm4
_mm_storeu_ps((pOut -= nHeight), row1);
000000013F40125C sub rdx,rbp
000000013F40125F shufps xmm5,xmm2,88h
000000013F401263 movups xmmword ptr [rdx],xmm0
_mm_storeu_ps((pOut -= nHeight), row0);
000000013F401266 sub rdx,rbp
000000013F401269 movups xmmword ptr [rdx],xmm5
000000013F40126C sub rax,1
000000013F401270 jne TransposeNonSquareMatrix_SSE+1A0h (013F401210h)
The compiler has split the loop variable into two, rax and r10. The final shuffles from the transpose macro are interleaved between the unaligned stores.
But often the matrix will be non-square, and so cannot be transposed in place. But let's look at transposing with a source and destination.
We will transpose in 4x4 blocks. Momentarily ignoring what is evenly divisible by 4:
for (size_t y = nHeight&~3; y >0;)
{
y -= 4;
for (size_t x = nWidth&~3; x>0;)
{
x -= 4;
__m128 row0;
__m128 row1;
__m128 row2;
__m128 row3;
row0 = _mm_loadu_ps(&pMatrixIn[(y + 0)*nWidth + x]);
row1 = _mm_loadu_ps(&pMatrixIn[(y + 1)*nWidth + x]);
row2 = _mm_loadu_ps(&pMatrixIn[(y + 2)*nWidth + x]);
row3 = _mm_loadu_ps(&pMatrixIn[(y + 3)*nWidth + x]);
_MM_TRANSPOSE4_PS(row0, row1, row2, row3);
_mm_storeu_ps(&pMatrixOut[(x + 0)*nHeight + y], row0);
_mm_storeu_ps(&pMatrixOut[(x + 1)*nHeight + y], row1);
_mm_storeu_ps(&pMatrixOut[(x + 2)*nHeight + y], row2);
_mm_storeu_ps(&pMatrixOut[(x + 3)*nHeight + y], row3);
}
}
This seems pretty tight...but that's a lot of multiplies. We can factor out most of the multiplies so that it's just pointer arithmetic and adds. On x64 the central loop can be accomplished without multiplies and without register spills.
Most of the matrix will be transposed by this set of loops. We will have to handle the rest in two other rectangular strips (making sure not to double-count the corner rectangle).
void TransposeNonSquareMatrix_SSE(const float *const pMatrixIn, float *const pMatrixOut, const size_t nWidth, const size_t nHeight)
{
for (size_t y = nHeight; y & 3;)
{
y--;
for (size_t x = nWidth; x-- > 0;)
{
pMatrixOut[x*nHeight + y] = pMatrixIn[y*nWidth + x];
}
}
if (nWidth & 3)
{
for (size_t y = nHeight&~3; y-- > 0;)
{
for (size_t x = nWidth; x & 3;)
{
x--;
pMatrixOut[x*nHeight + y] = pMatrixIn[y*nWidth + x];
}
}
}
const size_t nHeightMax = nHeight&~3;
const size_t nWidthMax = nWidth&~3;
for (size_t y = nHeightMax; y >0;)
{
y -= 4;
const float *const pIn0 = pMatrixIn + y*nWidth;
const float *const pIn1 = pIn0 + nWidth;
const float *const pIn2 = pIn1 + nWidth;
const float * const pIn3 = pIn2 + nWidth;
float * pOut = &pMatrixOut[nWidthMax*nHeight + y];
for (size_t x = nWidthMax; x>0;)
{
x -= 4;
__m128 row0;
__m128 row1;
__m128 row2;
__m128 row3;
row0 = _mm_loadu_ps(&pIn0[x]);
row1 = _mm_loadu_ps(&pIn1[x]);
row2 = _mm_loadu_ps(&pIn2[x]);
row3 = _mm_loadu_ps(&pIn3[x]);
_MM_TRANSPOSE4_PS(row0, row1, row2, row3);
_mm_storeu_ps((pOut -= nHeight), row3);
_mm_storeu_ps((pOut -= nHeight), row2);
_mm_storeu_ps((pOut -= nHeight), row1);
_mm_storeu_ps((pOut -= nHeight), row0);
}
}
}
The assembly compilation of the central loop is pretty tight:
{
x -= 4;
__m128 row0;
__m128 row1;
__m128 row2;
__m128 row3;
row0 = _mm_loadu_ps(&pIn0[x]);
row1 = _mm_loadu_ps(&pIn1[x]);
000000013F401210 movups xmm0,xmmword ptr [r10-10h]
row2 = _mm_loadu_ps(&pIn2[x]);
row3 = _mm_loadu_ps(&pIn3[x]);
_MM_TRANSPOSE4_PS(row0, row1, row2, row3);
_mm_storeu_ps((pOut -= nHeight), row3);
000000013F401215 sub rdx,rbp
000000013F401218 lea r10,[r10-10h]
000000013F40121C movups xmm4,xmmword ptr [r10+r15]
000000013F401221 movups xmm3,xmmword ptr [r10+r11]
000000013F401226 movups xmm1,xmmword ptr [r10+rcx]
000000013F40122B movaps xmm5,xmm4
000000013F40122E movaps xmm2,xmm3
000000013F401231 shufps xmm4,xmm0,0EEh
000000013F401235 shufps xmm5,xmm0,44h
row2 = _mm_loadu_ps(&pIn2[x]);
row3 = _mm_loadu_ps(&pIn3[x]);
_MM_TRANSPOSE4_PS(row0, row1, row2, row3);
_mm_storeu_ps((pOut -= nHeight), row3);
000000013F401239 movaps xmm0,xmm4
000000013F40123C shufps xmm3,xmm1,0EEh
000000013F401240 shufps xmm0,xmm3,0DDh
000000013F401244 movups xmmword ptr [rdx],xmm0
_mm_storeu_ps((pOut -= nHeight), row2);
000000013F401247 sub rdx,rbp
000000013F40124A shufps xmm2,xmm1,44h
000000013F40124E movaps xmm0,xmm5
000000013F401251 shufps xmm4,xmm3,88h
000000013F401255 shufps xmm0,xmm2,0DDh
000000013F401259 movups xmmword ptr [rdx],xmm4
_mm_storeu_ps((pOut -= nHeight), row1);
000000013F40125C sub rdx,rbp
000000013F40125F shufps xmm5,xmm2,88h
000000013F401263 movups xmmword ptr [rdx],xmm0
_mm_storeu_ps((pOut -= nHeight), row0);
000000013F401266 sub rdx,rbp
000000013F401269 movups xmmword ptr [rdx],xmm5
000000013F40126C sub rax,1
000000013F401270 jne TransposeNonSquareMatrix_SSE+1A0h (013F401210h)
The compiler has split the loop variable into two, rax and r10. The final shuffles from the transpose macro are interleaved between the unaligned stores.
Monday, January 2, 2017
Gauss-Jordan Elimination - Part2
Even after the last post, there's still a bunch of code necessary to get Gauss-Jordan elimination working. Most of it is pretty mundane nuts and bolts. But let's be thorough.
We need to multiply a vector in place by a runtime float:
void mul_ps_SSE2(float * const pA, size_t cAOrig, const float fac)
{
size_t uiA = cAOrig;
while (uiA & 3)
{
pA[--uiA] *= fac;
}
if (uiA)
{
__m128 mFac = _mm_set1_ps(fac);
do {
uiA -= 4;
_mm_storeu_ps(pA + uiA,_mm_mul_ps(_mm_loadu_ps(pA + uiA), mFac));
} while (uiA > 0);
}
}
We have a need to add a vector multiplied by a runtime float to another vector:
void addVectorTimesConstant_ps_SSE2(float * const pDest, const float * const pProd, size_t cAOrig, const float dum)
{
size_t uiA = cAOrig;
while (uiA & 3)
{
uiA--;
pDest[uiA] += pProd[uiA] * dum;
}
if (uiA)
{
const __m128 mDum = _mm_set1_ps(dum);
do {
uiA -= 4;
__m128 mDest = _mm_loadu_ps(pDest + uiA);
__m128 mProd = _mm_loadu_ps(pProd + uiA);
__m128 mResult = _mm_add_ps(mDest, _mm_mul_ps(mProd, mDum));
_mm_storeu_ps(pDest + uiA, mResult);
} while (uiA > 0);
}
}
While it's not huge, we can also SSE the assignment of a runtime int to a vector. Casting would allow this to be used for single precision floats.
void setVector_SSE2(int * const pDest, const size_t cAOrig, int iVal)
{
size_t uiA = cAOrig;
while (uiA & 3)
{
uiA--;
pDest[uiA] = iVal;
}
if (uiA)
{
const __m128i mVal= _mm_set1_epi32(iVal);
do {
uiA -= 4;
_mm_storeu_si128((__m128i *)pDest + uiA, mVal);
} while (uiA > 0);
}
}
Swapping the contents of two vectors isn't too difficult:
template <class T>
void swap(T &a, T&b)
{
T temp = a;
a = b;
b = temp;
}
void swapArray_SSE2(int * const pA, int * const pB, const size_t cAOrig)
{
size_t uiA = cAOrig;
if (uiA & 1)
{
uiA--;
swap(pA[uiA], pB[uiA]);
}
if (uiA & 2)
{
uiA-=2;
swap(* reinterpret_cast<__int64*>(pA+uiA), *reinterpret_cast<__int64*>(pA + uiA));
}
while (uiA > 0)
{
uiA -= 4;
__m128i mA = _mm_loadu_si128((__m128i*)(pA + uiA));
__m128i mB = _mm_loadu_si128((__m128i*)(pB + uiA));
_mm_storeu_si128((__m128i*)(pA + uiA), mB);
_mm_storeu_si128((__m128i*)(pB + uiA), mA);
}
}
void swapArray_SSE2(float * const pA, float * const pB, const size_t cAOrig)
{
swapArray_SSE2(reinterpret_cast<int*>(pA), reinterpret_cast<int*>(pB), cAOrig);
}
With these additional helper method it isn't too difficult to translate the implementation of gaussj from Numerical Recipes.
void gaussj_SSE2(float **a, const UINT n, float **b, const UINT m)
{
int *ipiv = (int*)malloc(n * sizeof(int));
UINT *indxc = (UINT*)malloc(n * sizeof(UINT));
UINT *indxr = (UINT*)malloc(n * sizeof(UINT));
setVector_SSE2(ipiv, n, 0x7fffffff);
for (size_t i = n; i-->0; i++)
{
float big = 0;
size_t irow = 0, icol = 0;
for (size_t j = n; j-- >0; )
{
if (ipiv[j])
{
const float *pfBestInRow = findLargestMaskedFloat_SSE2(ipiv, a[j], (UINT)n);
float fBestValInRow = fabsf(*pfBestInRow);
if (fBestValInRow > big)
{
big = fBestValInRow;
irow = j;
icol = (size_t)(pfBestInRow - a[j]);
}
}
}
ipiv[icol] = 0;
if (irow != icol)
{
swapArray_SSE2(a[irow], a[icol], n);
swapArray_SSE2(b[irow], b[icol], m);
}
indxr[i] = (UINT)irow;
indxc[i] = (UINT)icol;
if (0 == a[icol][icol])
{
assert(0 && "gaussj : Singular Matrix");
return;
}
float pivinv = 1.0f / a[icol][icol];
a[icol][icol] = 1;
mul_ps_SSE2(a[icol], n, pivinv);
mul_ps_SSE2(b[icol], m, pivinv);
for (size_t ll = n; ll-- >0; )
{
if (ll != icol)
{
float dum = -a[ll][icol];
a[ll][icol] = 0;
addVectorTimesConstant_ps_SSE2(a[ll], a[icol], n, dum);
addVectorTimesConstant_ps_SSE2(b[ll], b[icol], m, dum);
}
}
}
for (size_t l = n ; l-- > 0; )
{
if (indxr[l] != indxc[l])
{
for (size_t k = 0; k < n; k++)
{
swap(a[k][indxr[l]], a[k][indxc[l]]);
}
}
}
free(indxc);
free(indxr);
free(ipiv);
}
We need to multiply a vector in place by a runtime float:
void mul_ps_SSE2(float * const pA, size_t cAOrig, const float fac)
{
size_t uiA = cAOrig;
while (uiA & 3)
{
pA[--uiA] *= fac;
}
if (uiA)
{
__m128 mFac = _mm_set1_ps(fac);
do {
uiA -= 4;
_mm_storeu_ps(pA + uiA,_mm_mul_ps(_mm_loadu_ps(pA + uiA), mFac));
} while (uiA > 0);
}
}
We have a need to add a vector multiplied by a runtime float to another vector:
void addVectorTimesConstant_ps_SSE2(float * const pDest, const float * const pProd, size_t cAOrig, const float dum)
{
size_t uiA = cAOrig;
while (uiA & 3)
{
uiA--;
pDest[uiA] += pProd[uiA] * dum;
}
if (uiA)
{
const __m128 mDum = _mm_set1_ps(dum);
do {
uiA -= 4;
__m128 mDest = _mm_loadu_ps(pDest + uiA);
__m128 mProd = _mm_loadu_ps(pProd + uiA);
__m128 mResult = _mm_add_ps(mDest, _mm_mul_ps(mProd, mDum));
_mm_storeu_ps(pDest + uiA, mResult);
} while (uiA > 0);
}
}
While it's not huge, we can also SSE the assignment of a runtime int to a vector. Casting would allow this to be used for single precision floats.
void setVector_SSE2(int * const pDest, const size_t cAOrig, int iVal)
{
size_t uiA = cAOrig;
while (uiA & 3)
{
uiA--;
pDest[uiA] = iVal;
}
if (uiA)
{
const __m128i mVal= _mm_set1_epi32(iVal);
do {
uiA -= 4;
_mm_storeu_si128((__m128i *)pDest + uiA, mVal);
} while (uiA > 0);
}
}
Swapping the contents of two vectors isn't too difficult:
template <class T>
void swap(T &a, T&b)
{
T temp = a;
a = b;
b = temp;
}
void swapArray_SSE2(int * const pA, int * const pB, const size_t cAOrig)
{
size_t uiA = cAOrig;
if (uiA & 1)
{
uiA--;
swap(pA[uiA], pB[uiA]);
}
if (uiA & 2)
{
uiA-=2;
swap(* reinterpret_cast<__int64*>(pA+uiA), *reinterpret_cast<__int64*>(pA + uiA));
}
while (uiA > 0)
{
uiA -= 4;
__m128i mA = _mm_loadu_si128((__m128i*)(pA + uiA));
__m128i mB = _mm_loadu_si128((__m128i*)(pB + uiA));
_mm_storeu_si128((__m128i*)(pA + uiA), mB);
_mm_storeu_si128((__m128i*)(pB + uiA), mA);
}
}
void swapArray_SSE2(float * const pA, float * const pB, const size_t cAOrig)
{
swapArray_SSE2(reinterpret_cast<int*>(pA), reinterpret_cast<int*>(pB), cAOrig);
}
With these additional helper method it isn't too difficult to translate the implementation of gaussj from Numerical Recipes.
void gaussj_SSE2(float **a, const UINT n, float **b, const UINT m)
{
int *ipiv = (int*)malloc(n * sizeof(int));
UINT *indxc = (UINT*)malloc(n * sizeof(UINT));
UINT *indxr = (UINT*)malloc(n * sizeof(UINT));
setVector_SSE2(ipiv, n, 0x7fffffff);
for (size_t i = n; i-->0; i++)
{
float big = 0;
size_t irow = 0, icol = 0;
for (size_t j = n; j-- >0; )
{
if (ipiv[j])
{
const float *pfBestInRow = findLargestMaskedFloat_SSE2(ipiv, a[j], (UINT)n);
float fBestValInRow = fabsf(*pfBestInRow);
if (fBestValInRow > big)
{
big = fBestValInRow;
irow = j;
icol = (size_t)(pfBestInRow - a[j]);
}
}
}
ipiv[icol] = 0;
if (irow != icol)
{
swapArray_SSE2(a[irow], a[icol], n);
swapArray_SSE2(b[irow], b[icol], m);
}
indxr[i] = (UINT)irow;
indxc[i] = (UINT)icol;
if (0 == a[icol][icol])
{
assert(0 && "gaussj : Singular Matrix");
return;
}
float pivinv = 1.0f / a[icol][icol];
a[icol][icol] = 1;
mul_ps_SSE2(a[icol], n, pivinv);
mul_ps_SSE2(b[icol], m, pivinv);
for (size_t ll = n; ll-- >0; )
{
if (ll != icol)
{
float dum = -a[ll][icol];
a[ll][icol] = 0;
addVectorTimesConstant_ps_SSE2(a[ll], a[icol], n, dum);
addVectorTimesConstant_ps_SSE2(b[ll], b[icol], m, dum);
}
}
}
for (size_t l = n ; l-- > 0; )
{
if (indxr[l] != indxc[l])
{
for (size_t k = 0; k < n; k++)
{
swap(a[k][indxr[l]], a[k][indxc[l]]);
}
}
}
free(indxc);
free(indxr);
free(ipiv);
}
Subscribe to:
Posts (Atom)