diff --git a/GPU/GPUTracking/Definitions/GPUSettingsList.h b/GPU/GPUTracking/Definitions/GPUSettingsList.h index e3327f2f8b661..0dc293ae7be38 100644 --- a/GPU/GPUTracking/Definitions/GPUSettingsList.h +++ b/GPU/GPUTracking/Definitions/GPUSettingsList.h @@ -113,7 +113,7 @@ AddOptionRTC(maxTimeBinAboveThresholdIn1000Bin, uint16_t, 500, "", 0, "Except pa AddOptionRTC(maxConsecTimeBinAboveThreshold, uint16_t, 200, "", 0, "Except pad from cluster finding if number of consecutive charges in a fragment is above this baseline (disable = 0)") AddOptionRTC(noisyPadSaturationThreshold, uint16_t, 700, "", 0, "Threshold where a timebin is considered saturated, disabling the noisy pad check for that pad") AddOptionRTC(hipTailFilter, uint8_t, 0, "", 0, "Enable Highly Ionising Particle tail filter in CheckPadBaseline (0 = disable, 1 = filter tails)") -AddOptionRTC(hipTailFilterMinimum, uint16_t, 1024, "", 0, "Thread signal above this minimum as saturated") +AddOptionRTC(hipTailFilterMinimum, uint16_t, 1023, "", 0, "Thread signal above this minimum as saturated") AddOptionRTC(hipTailFilterThreshold, uint16_t, 100, "", 0, "Threshold that must be exceeded for a timebin to be counted towards Highly Ionising Particle tail") AddOptionRTC(hipTailFilterAlpha, float, 0.5f, "", 0, "Smoothing factor for the exponential Highly Ionising Particle tail filter") AddOptionRTC(occupancyMapTimeBins, uint16_t, 16, "", 0, "Number of timebins per histogram bin of occupancy map (0 = disable occupancy map)") diff --git a/GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx b/GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx index 76dcd54a89289..f86c4ba18cee0 100644 --- a/GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx +++ b/GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx @@ -1144,16 +1144,12 @@ int32_t GPUChainTracking::RunTPCClusterizer(bool synchronizeOutput) // TODO Add some warning when re enabling pad filter with this flag, so it's not just silently enabled when disabling was requested checkForNoisyPads |= rec()->GetParam().rec.tpc.hipTailFilter; - if (rec()->GetParam().rec.tpc.hipTailFilter && !doGPU) { - GPUError("HIP tail filter enabled, but this is currently not supported on CPU"); - } - if (checkForNoisyPads) { if (rec()->GetParam().rec.tpc.hipTailFilter) { runKernel({GetGridAutoStep(lane, RecoStep::TPCClusterFinding)}, clustererShadow.mPhipTailsByRow, GPUTPCGeometry::NROWS * sizeof(*clustererShadow.mPhipTailsByRow) * GPUTPCCFHIPClusterizer::MaxHIPTailsPerRow); runKernel({GetGridAutoStep(lane, RecoStep::TPCClusterFinding)}, clustererShadow.mPnHIPTails, GPUTPCGeometry::NROWS * sizeof(*clustererShadow.mPnHIPTails)); } - const int32_t nBlocks = GPUTPCCFCheckPadBaseline::GetNBlocks(doGPU); + const int32_t nBlocks = GPUTPCGeometry::NROWS; runKernel({GetGridBlk(nBlocks, lane), {iSector}}); getKernelTimer(RecoStep::TPCClusterFinding, iSector, TPC_REAL_PADS_IN_SECTOR * fragment.lengthWithoutOverlap() * sizeof(PackedCharge), false); diff --git a/GPU/GPUTracking/TPCClusterFinder/GPUTPCCFCheckPadBaseline.cxx b/GPU/GPUTracking/TPCClusterFinder/GPUTPCCFCheckPadBaseline.cxx index 5a5fb478fe065..7470f86877490 100644 --- a/GPU/GPUTracking/TPCClusterFinder/GPUTPCCFCheckPadBaseline.cxx +++ b/GPU/GPUTracking/TPCClusterFinder/GPUTPCCFCheckPadBaseline.cxx @@ -22,6 +22,7 @@ #ifndef GPUCA_GPUCODE #include "MCLabelAccumulator.h" #include "utils/VcShim.h" +#include #endif #if 0 @@ -344,7 +345,7 @@ GPUd() void GPUTPCCFCheckPadBaseline::CheckBaselineGPU(int32_t nBlocks, int32_t GPUbarrier(); - } // if (hipTriggerFound) + } // if (hasHIPTrigger) } // for (uint16_t t = firstTB; t < lastTB; t += NumOfCachedTBs) @@ -374,60 +375,264 @@ GPUd() void GPUTPCCFCheckPadBaseline::CheckBaselineGPU(int32_t nBlocks, int32_t GPUd() void GPUTPCCFCheckPadBaseline::CheckBaselineCPU(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUSharedMemory& smem, processorType& clusterer) { #ifndef GPUCA_GPUCODE - const CfFragment& fragment = clusterer.mPmemory->fragment; - CfArray2D chargeMap(reinterpret_cast(clusterer.mPchargeMap)); - - CfChargePos basePos(iBlock * PadsPerCacheline, 0); - - constexpr GPUTPCGeometry geo; - if (basePos.pad() >= geo.NPads(basePos.row())) { + if (iBlock >= (int32_t)GPUTPCGeometry::NROWS) { return; } - constexpr size_t ElemsInTileRow = (size_t)TilingLayout>::WidthInTiles * TimebinsPerCacheline * PadsPerCacheline; + constexpr GPUTPCGeometry geo; + const int32_t row = iBlock; + const int32_t nPads = geo.NPads(row); + const int32_t nVecPads = (nPads + PadsPerCacheline - 1) / PadsPerCacheline; + + const CfFragment& fragment = clusterer.mPmemory->fragment; + const bool hipFilterOn = clusterer.Param().rec.tpc.hipTailFilter; + const Charge hipTailThreshold = clusterer.Param().rec.tpc.hipTailFilterThreshold; + const Charge hipTailFilterAlpha = clusterer.Param().rec.tpc.hipTailFilterAlpha; + auto* nHIPTails = clusterer.mPnHIPTails; + auto* hipTails = GetHIPTails(clusterer, row); + + CfArray2D chargeMap(reinterpret_cast(clusterer.mPchargeMap)); using UShort8 = Vc::fixed_size_simd; + using Short8 = Vc::fixed_size_simd; using Charge8 = Vc::fixed_size_simd; - UShort8 totalCharges{Vc::Zero}; - UShort8 consecCharges{Vc::Zero}; - UShort8 maxConsecCharges{Vc::Zero}; - Charge8 maxCharge{Vc::Zero}; - - tpccf::TPCFragmentTime t = fragment.firstNonOverlapTimeBin(); - - // Access packed charges as raw integers. We throw away the PackedCharge type here to simplify vectorization. - const uint16_t* packedChargeStart = reinterpret_cast(&chargeMap[basePos.delta({0, t})]); - - for (; t < fragment.lastNonOverlapTimeBin(); t += TimebinsPerCacheline) { - for (tpccf::TPCFragmentTime localtime = 0; localtime < TimebinsPerCacheline; localtime++) { - const UShort8 packedCharges{packedChargeStart + PadsPerCacheline * localtime, Vc::Aligned}; - const UShort8::mask_type isCharge = packedCharges != 0; - - if (isCharge.isNotEmpty()) { - totalCharges(isCharge)++; - consecCharges += 1; - consecCharges(not isCharge) = 0; - maxConsecCharges = Vc::max(consecCharges, maxConsecCharges); - - // Manually unpack charges to float. - // Duplicated from PackedCharge::unpack to generate vectorized code: - // Charge unpack() const { return Charge(mVal & ChargeMask) / Charge(1 << DecimalBits); } - // Note that PackedCharge has to cut off the highest 2 bits via ChargeMask as they are used for flags by the cluster finder - // and are not part of the charge value. We can skip this step because the cluster finder hasn't run yet - // and thus these bits are guarenteed to be zero. - const Charge8 unpackedCharges = Charge8(packedCharges) / Charge(1 << PackedCharge::DecimalBits); - maxCharge = Vc::max(maxCharge, unpackedCharges); - } else { - consecCharges = 0; - } + std::vector totalChargesV(nVecPads, UShort8{Vc::Zero}); + std::vector consecChargesV(nVecPads, UShort8{Vc::Zero}); + std::vector maxConsecChargesV(nVecPads, UShort8{Vc::Zero}); + std::vector maxChargeV(nVecPads, Charge8{Vc::Zero}); + + std::vector localHipTbV(nVecPads, -1); + std::vector broadcastHipTbV(nVecPads, -1); + std::vector aboveThresholdStartV(nVecPads, -1); + std::vector activeHIPTailStartV(nVecPads, -1); + std::vector activeHIPTailEndV(nVecPads, -1); + std::vector tailFilterChargeV(nVecPads, Charge8{Vc::Zero}); + + for (int16_t t = 0; t < fragment.length; t += NumOfCachedTBs) { + + bool hasAnyTrigger = false; + + // Run actual noisy pad filter and look for HIP trigger + for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) { + + auto totalCharges = totalChargesV[iVecPad]; + auto consecCharges = consecChargesV[iVecPad]; + auto maxConsecCharges = maxConsecChargesV[iVecPad]; + auto maxCharge = maxChargeV[iVecPad]; + + auto hipTb = Short8(-1); + auto aboveThresholdStart = aboveThresholdStartV[iVecPad]; + auto activeHIPTailStart = activeHIPTailStartV[iVecPad]; + auto activeHIPTailEnd = activeHIPTailEndV[iVecPad]; + auto tailFilterCharge = tailFilterChargeV[iVecPad]; + + const CfChargePos basePos(row, iVecPad * PadsPerCacheline, t); + + for (tpccf::TPCFragmentTime localtime = 0; localtime < NumOfCachedTBs; localtime++) { + + const uint16_t* packedChargeStart = reinterpret_cast(&chargeMap[basePos.delta({0, localtime})]); + const UShort8 packedCharges = t + localtime < fragment.length + ? UShort8{packedChargeStart, Vc::Aligned} + : UShort8{Vc::Zero}; + const auto isCharge = packedCharges != 0; + + const auto unpackedCharges = Charge8(packedCharges) / Charge(1 << PackedCharge::DecimalBits); + + if (isCharge.isNotEmpty()) { + totalCharges(isCharge)++; + consecCharges += 1; + consecCharges(not isCharge) = 0; + maxConsecCharges = Vc::max(consecCharges, maxConsecCharges); + + // Manually unpack charges to float. + // Duplicated from PackedCharge::unpack to generate vectorized code: + // Charge unpack() const { return Charge(mVal & ChargeMask) / Charge(1 << DecimalBits); } + // Note that PackedCharge has to cut off the highest 2 bits via ChargeMask as they are used for flags by the cluster finder + // and are not part of the charge value. We can skip this step because the cluster finder hasn't run yet + // and thus these bits are guarenteed to be zero. + maxCharge = Vc::max(maxCharge, unpackedCharges); + + const auto aboveRisingEdge = unpackedCharges >= hipTailThreshold; + const auto startRisingEdge = aboveRisingEdge && aboveThresholdStart < 0; + aboveThresholdStart(startRisingEdge) = t + localtime; + aboveThresholdStart(!aboveRisingEdge) = -1; + + const auto hasNewTrigger = hipTb < 0 && unpackedCharges >= Charge(MaxADC); + hipTb(hasNewTrigger) = aboveThresholdStart; + hasAnyTrigger |= hasNewTrigger.isNotEmpty(); + } else { + consecCharges = 0; + aboveThresholdStart = -1; + } + + const auto tailOpen = activeHIPTailStart > -1 && activeHIPTailEnd < 0; + tailFilterCharge(tailOpen) = tailFilterCharge + hipTailFilterAlpha * (unpackedCharges - tailFilterCharge); + activeHIPTailEnd(tailOpen && tailFilterCharge < hipTailThreshold) = t + localtime; + } // for (tpccf::TPCFragmentTime localtime = 0; localtime < TimebinsPerCacheline; localtime++) + + totalChargesV[iVecPad] = totalCharges; + consecChargesV[iVecPad] = consecCharges; + maxConsecChargesV[iVecPad] = maxConsecCharges; + maxChargeV[iVecPad] = maxCharge; + + localHipTbV[iVecPad] = hipTb; + aboveThresholdStartV[iVecPad] = aboveThresholdStart; + activeHIPTailStartV[iVecPad] = activeHIPTailStart; + activeHIPTailEndV[iVecPad] = activeHIPTailEnd; + tailFilterChargeV[iVecPad] = tailFilterCharge; + + } // for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) + + if (hasAnyTrigger) { + broadcastHipTbV = localHipTbV; } - packedChargeStart += ElemsInTileRow; - } + // Broadcast trigger times to neighboring pads across the whole + for (int16_t iVecPad = 0; iVecPad < nVecPads && hasAnyTrigger; iVecPad++) { + + const auto hipTb = localHipTbV[iVecPad]; + + const auto hasHipTrigger = hipTb > -1; + if (hasHipTrigger.isNotEmpty()) [[unlikely]] { + + // TODO: This could be vectorised, but doesn't seem necessary + for (uint16_t p = 0; p < PadsPerCacheline; p++) { + if (hasHipTrigger[p]) { + const int16_t pad = iVecPad * PadsPerCacheline + p; + const int16_t neighborSt = CAMath::Max(0, pad - SSClusterPadWidth); + const int16_t neighborEnd = CAMath::Min(nPads, pad + SSClusterPadWidth + 1); + for (int16_t np = neighborSt; np < neighborEnd; np++) { + if (np == pad) { + continue; + } + const auto pv = np / PadsPerCacheline; + const auto pi = np % PadsPerCacheline; + // GPU keeps a pad's own trigger time; only pads without a local trigger inherit from neighbors. + if (localHipTbV[pv][pi] < 0) { + broadcastHipTbV[pv][pi] = CAMath::Max(hipTb[p], broadcastHipTbV[pv][pi]); + } + } // for (int16_t np = neighborSt; np < neighborEnd; np++) + } // if (hasHipTrigger[p]) { + } // for (uint16_t p = 0; p < PadsPerCacheline; p++) + } // if (hasHipTrigger.isNotEmpty()) + } // for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) + + // Close old tails for all pads, open new tails in case of overlap + for (int16_t iVecPad = 0; iVecPad < nVecPads && hasAnyTrigger; iVecPad++) { + + auto hipTb = broadcastHipTbV[iVecPad]; + auto aboveThresholdStart = aboveThresholdStartV[iVecPad]; + auto activeHIPTailStart = activeHIPTailStartV[iVecPad]; + auto activeHIPTailEnd = activeHIPTailEndV[iVecPad]; + auto tailFilterCharge = tailFilterChargeV[iVecPad]; + + const auto shouldCloseTail = hipTb > -1 && activeHIPTailStart > -1; + activeHIPTailEnd(shouldCloseTail && activeHIPTailEnd < 0) = hipTb; + + // Closing tails will store them to global memory and zero the range + // So it's enough to disable this part to fully disable the tail filter + if (hipFilterOn && shouldCloseTail.isNotEmpty()) { + for (int16_t p = 0; p < PadsPerCacheline; p++) { + const int16_t pad = iVecPad * PadsPerCacheline + p; + if (shouldCloseTail[p] && pad < nPads) { + Charge tailQtot = 0; + Charge tailQMax = 0; + for (int16_t tt = activeHIPTailStart[p]; tt < activeHIPTailEnd[p]; tt++) { + const CfChargePos basePos(row, iVecPad * PadsPerCacheline, 0); + const auto pos = basePos.delta({p, tt}); + const auto q = chargeMap[pos].unpack(); + tailQtot += q; + tailQMax = CAMath::Max(tailQMax, q); + chargeMap[pos] = PackedCharge{0}; + } + + if (activeHIPTailEnd[p] > activeHIPTailStart[p]) { // Prune empty tails + const auto tailIdx = CAMath::AtomicAdd(&nHIPTails[row], 1) + 1; + if (tailIdx < GPUTPCCFHIPTailConnector::MaxHIPTailsPerRow) { + hipTails[tailIdx] = { + .iPrev = 0, + .iNext = 0, + .pad = uint16_t(pad), + .tailStart = uint16_t(activeHIPTailStart[p]), + .tailEnd = uint16_t(activeHIPTailEnd[p]), + .qTot = tailQtot, + .qMax = tailQMax, + }; + } + } + + } // if (shouldCloseTail[p] && pad < nPads) + } // for (uint16_t p = 0; p < PadsPerCacheline; p++) + } // if (shouldCloseThipFilterOn && shouldCloseTail.isNotEmpty()) + + activeHIPTailStart(hipTb > -1) = hipTb; + activeHIPTailEnd(hipTb > -1) = -1; + tailFilterCharge(hipTb > -1) = MaxADC; + + aboveThresholdStartV[iVecPad] = aboveThresholdStart; + activeHIPTailStartV[iVecPad] = activeHIPTailStart; + activeHIPTailEndV[iVecPad] = activeHIPTailEnd; + tailFilterChargeV[iVecPad] = tailFilterCharge; + + } // for (int32_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) + } // for (auto t = 0; t < fragment.length; t += TimebinsPerCacheline) + + // Close old tails for all pads, open new tails in case of overlap + for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) { + + auto activeHIPTailStart = activeHIPTailStartV[iVecPad]; + auto activeHIPTailEnd = activeHIPTailEndV[iVecPad]; + + const auto shouldCloseTail = activeHIPTailStart > -1; + activeHIPTailEnd(shouldCloseTail && activeHIPTailEnd < 0) = fragment.length; + + if (hipFilterOn && shouldCloseTail.isNotEmpty()) { + for (int16_t p = 0; p < PadsPerCacheline; p++) { + const int16_t pad = iVecPad * PadsPerCacheline + p; + if (shouldCloseTail[p] && pad < nPads) { + Charge tailQtot = 0; + Charge tailQMax = 0; + for (int16_t tt = activeHIPTailStart[p]; tt < activeHIPTailEnd[p]; tt++) { + const CfChargePos basePos(row, iVecPad * PadsPerCacheline, 0); + const auto pos = basePos.delta({p, tt}); + const auto q = chargeMap[pos].unpack(); + tailQtot += q; + tailQMax = CAMath::Max(tailQMax, q); + chargeMap[pos] = PackedCharge{0}; + } + + if (activeHIPTailEnd[p] > activeHIPTailStart[p]) { // Prune empty tails + const auto tailIdx = CAMath::AtomicAdd(&nHIPTails[row], 1) + 1; + if (tailIdx < GPUTPCCFHIPTailConnector::MaxHIPTailsPerRow) { + hipTails[tailIdx] = { + .iPrev = 0, + .iNext = 0, + .pad = uint16_t(pad), + .tailStart = uint16_t(activeHIPTailStart[p]), + .tailEnd = uint16_t(activeHIPTailEnd[p]), + .qTot = tailQtot, + .qMax = tailQMax, + }; + } + } + + } // if (shouldCloseTail[p] && pad < nPads) + } // for (uint16_t p = 0; p < PadsPerCacheline; p++) + } // if (hipFilterOn && shouldCloseTail.isNotEmpty()) + } // for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) + + for (int32_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) { + + const UShort8 totalCharges = totalChargesV[iVecPad]; + const UShort8 maxConsecCharges = maxConsecChargesV[iVecPad]; + const Charge8 maxCharge = maxChargeV[iVecPad]; + + const CfChargePos basePos(row, iVecPad * PadsPerCacheline, 0); - for (tpccf::Pad localpad = 0; localpad < PadsPerCacheline; localpad++) { - updatePadBaseline(basePos.gpad + localpad, clusterer, totalCharges[localpad], maxConsecCharges[localpad], maxCharge[localpad]); + for (tpccf::Pad localpad = 0; localpad < PadsPerCacheline; localpad++) { + updatePadBaseline(basePos.gpad + localpad, clusterer, totalCharges[localpad], maxConsecCharges[localpad], maxCharge[localpad]); + } } #endif } @@ -463,16 +668,23 @@ GPUd() void GPUTPCCFHIPTailConnector::Thread<0>(int32_t nBlocks, int32_t nThread #ifdef GPUCA_DETERMINISTIC_MODE // Races in tail comparisons and atomic swap can lead to slightly different clusters. // So need a sequential fallback for deterministic mode - if (iThread > 0) { - return; - } - nThreads = 1; GPUCommonAlgorithm::sortInBlock(tails + 1, tails + nTails + 1, [](auto&& t1, auto&& t2) { if (t1.pad != t2.pad) { return t1.pad < t2.pad; + } else if (t1.tailStart != t2.tailStart) { + return t1.tailStart < t2.tailStart; + } else if (t1.tailEnd != t2.tailEnd) { + return t1.tailEnd < t2.tailEnd; + } else if (t1.qTot != t2.qTot) { + return t1.qTot < t2.qTot; + } else { + return t1.qMax < t2.qMax; } - return t1.tailStart < t2.tailStart; }); + if (iThread > 0) { + return; + } + nThreads = 1; #endif for (uint32_t iTail = iThread + 1; iTail <= nTails; iTail += nThreads) { diff --git a/GPU/GPUTracking/TPCClusterFinder/GPUTPCCFCheckPadBaseline.h b/GPU/GPUTracking/TPCClusterFinder/GPUTPCCFCheckPadBaseline.h index c2c5a1e339256..08e6110ca2373 100644 --- a/GPU/GPUTracking/TPCClusterFinder/GPUTPCCFCheckPadBaseline.h +++ b/GPU/GPUTracking/TPCClusterFinder/GPUTPCCFCheckPadBaseline.h @@ -126,14 +126,6 @@ class GPUTPCCFCheckPadBaseline : public GPUKernelTemplate return gpudatatypes::RecoStep::TPCClusterFinding; } - static int32_t GetNBlocks(bool isGPU) - { - // Important to exclude rightmost padding from Pad Filter. - // There's nothing to filter there and padding is counted as start of a row, so it causes an overflow in the row count. - const int32_t nBlocksCPU = (TPC_CLUSTERER_STRIDED_PAD_COUNT - GPUCF_PADDING_PAD) / PadsPerCacheline; - return isGPU ? GPUTPCGeometry::NROWS : nBlocksCPU; - } - template GPUd() static void Thread(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUSharedMemory& smem, processorType& clusterer); diff --git a/GPU/GPUTracking/utils/VcShim.h b/GPU/GPUTracking/utils/VcShim.h index 2bbc1d471bbbb..b51210100c5a6 100644 --- a/GPU/GPUTracking/utils/VcShim.h +++ b/GPU/GPUTracking/utils/VcShim.h @@ -27,6 +27,7 @@ #include #include #include +#include namespace Vc { @@ -59,6 +60,7 @@ class WriteMaskVector V& mVec; public: + using vector_type = V; using value_type = typename V::value_type; WriteMaskVector(V& v, const M& m) : mMask(m), mVec(v) {} @@ -78,6 +80,15 @@ class WriteMaskVector } return *this; } + + WriteMaskVector& operator=(const vector_type& v) + { + for (size_t i = 0; i < mVec.size(); i++) { + if (mMask[i]) + mVec[i] = v[i]; + } + return *this; + } }; inline void prefetchMid(const void*) {} @@ -86,7 +97,7 @@ inline void prefetchForOneRead(const void*) {} } // namespace Common -template +template class fixed_size_simd_mask { private: @@ -95,7 +106,7 @@ class fixed_size_simd_mask public: bool isNotEmpty() const { return mData.any(); } - std::bitset::reference operator[](size_t i) { return mData[i]; } + typename std::bitset::reference operator[](size_t i) { return mData[i]; } bool operator[](size_t i) const { return mData[i]; } fixed_size_simd_mask operator!() const @@ -104,6 +115,13 @@ class fixed_size_simd_mask o.mData.flip(); return o; } + + fixed_size_simd_mask operator&&(const fixed_size_simd_mask& o) const + { + auto r = *this; + r.mData &= o.mData; + return r; + } }; template @@ -115,7 +133,7 @@ class fixed_size_simd public: using vector_type = std::array; using value_type = T; - using mask_type = fixed_size_simd_mask; + using mask_type = fixed_size_simd_mask; static constexpr size_t size() { return N; } @@ -130,6 +148,8 @@ class fixed_size_simd fixed_size_simd(const T* d, AlignedTag) { std::copy_n(d, N, mData.begin()); } + fixed_size_simd(const T& x) { mData.fill(x); } + T& operator[](size_t i) { return mData[i]; } const T& operator[](size_t i) const { return mData[i]; } @@ -149,6 +169,44 @@ class fixed_size_simd return *this; } + template + fixed_size_simd& operator+=(const fixed_size_simd& v) + { + for (size_t i = 0; i < N; i++) + mData[i] += v[i]; + return *this; + } + + fixed_size_simd& operator-=(const T& v) + { + for (auto& x : mData) + x -= v; + return *this; + } + + template + fixed_size_simd& operator-=(const fixed_size_simd& v) + { + for (size_t i = 0; i < N; i++) + mData[i] -= v[i]; + return *this; + } + + fixed_size_simd& operator*=(const T& v) + { + for (auto& x : mData) + x *= v; + return *this; + } + + template + fixed_size_simd& operator*=(const fixed_size_simd& v) + { + for (size_t i = 0; i < N; i++) + mData[i] *= v[i]; + return *this; + } + fixed_size_simd& operator/=(const T& v) { for (auto& x : mData) @@ -156,10 +214,12 @@ class fixed_size_simd return *this; } - fixed_size_simd operator/(const T& v) const + template + fixed_size_simd& operator/=(const fixed_size_simd& v) { - auto x = *this; - return x /= v; + for (size_t i = 0; i < N; i++) + mData[i] /= v[i]; + return *this; } mask_type operator==(const T& v) const @@ -172,10 +232,124 @@ class fixed_size_simd mask_type operator!=(const T& v) const { return !(*this == v); } + mask_type operator>(const T& v) const + { + mask_type m; + for (size_t i = 0; i < N; i++) + m[i] = mData[i] > v; + return m; + } + + mask_type operator>=(const T& v) const + { + mask_type m; + for (size_t i = 0; i < N; i++) + m[i] = mData[i] >= v; + return m; + } + + mask_type operator<(const T& v) const + { + mask_type m; + for (size_t i = 0; i < N; i++) + m[i] = mData[i] < v; + return m; + } + friend vector_type& internal_data<>(fixed_size_simd& x); friend const vector_type& internal_data<>(const fixed_size_simd& x); }; +template +struct is_fixed_size_simd : std::false_type { +}; + +template +struct is_fixed_size_simd> : std::true_type { +}; + +template +using EnableIfScalar = typename std::enable_if_t< + !is_fixed_size_simd>::value && std::is_convertible_v, int>; + +template +fixed_size_simd operator+(fixed_size_simd a, const fixed_size_simd& b) +{ + return a += b; +} + +template = 0> +fixed_size_simd operator+(fixed_size_simd a, const S& b) +{ + return a += static_cast(b); +} + +template = 0> +fixed_size_simd operator+(const S& a, fixed_size_simd b) +{ + return b += static_cast(a); +} + +template +fixed_size_simd operator-(fixed_size_simd a, const fixed_size_simd& b) +{ + return a -= b; +} + +template = 0> +fixed_size_simd operator-(fixed_size_simd a, const S& b) +{ + return a -= static_cast(b); +} + +template = 0> +fixed_size_simd operator-(const S& a, const fixed_size_simd& b) +{ + fixed_size_simd o; + for (size_t i = 0; i < N; i++) + o[i] = static_cast(a) - b[i]; + return o; +} + +template +fixed_size_simd operator*(fixed_size_simd a, const fixed_size_simd& b) +{ + return a *= b; +} + +template = 0> +fixed_size_simd operator*(fixed_size_simd a, const S& b) +{ + return a *= static_cast(b); +} + +template = 0> +fixed_size_simd operator*(const S& a, fixed_size_simd b) +{ + return b *= static_cast(a); +} + +template +fixed_size_simd operator/(fixed_size_simd a, const fixed_size_simd& b) +{ + return a /= b; +} + +template = 0> +fixed_size_simd operator/(fixed_size_simd a, const S& b) +{ + return a /= static_cast(b); +} + +template = 0> +fixed_size_simd operator/(const S& a, const fixed_size_simd& b) +{ + fixed_size_simd o; + for (size_t i = 0; i < N; i++) + o[i] = static_cast(a) / b[i]; + return o; +} + template V max(const V& a, const V& b) {