From 6717005529e8f8107ccf0eabe886e75031ae6f2c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 25 Sep 2025 03:16:49 +0200 Subject: [PATCH] Chunk GPU divisor cycle generation to avoid CPU fallback --- PerfectNumbers.Core/DivisorCycleCache.cs | 121 +++++++++++++++++++++-- 1 file changed, 111 insertions(+), 10 deletions(-) diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 34cac64..c7378fa 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -69,10 +69,13 @@ public sealed class DivisorCycleCache } private const int CycleGenerationBatchSize = 262_144; + private const int GpuCycleStepsPerInvocation = 524_288; + private const byte ByteZero = 0; + private const byte ByteOne = 1; private readonly object _sync = new(); private readonly ConcurrentDictionary> _pending = new(); - private readonly ConcurrentDictionary, ArrayView1D>> _gpuKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>> _gpuKernelCache = new(AcceleratorReferenceComparer.Instance); private CycleBlock _baseBlock = null!; private CycleBlock? _activeBlock; private CycleBlock? _prefetchedBlock; @@ -328,42 +331,102 @@ public sealed class DivisorCycleCache try { Accelerator accelerator = gpuLease.Accelerator; - Action, ArrayView1D> kernel = _gpuKernelCache.GetOrAdd(accelerator, LoadKernel); + Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> kernel = _gpuKernelCache.GetOrAdd(accelerator, LoadKernel); int length = destination.Length; int chunkCapacity = Math.Min(length, CycleGenerationBatchSize); ulong[] divisors = ArrayPool.Shared.Rent(chunkCapacity); + ulong[] powScratch = ArrayPool.Shared.Rent(chunkCapacity); + ulong[] orderScratch = ArrayPool.Shared.Rent(chunkCapacity); + byte[] statusScratch = ArrayPool.Shared.Rent(chunkCapacity); try { using MemoryBuffer1D divisorBuffer = accelerator.Allocate1D(chunkCapacity); + using MemoryBuffer1D powBuffer = accelerator.Allocate1D(chunkCapacity); + using MemoryBuffer1D orderBuffer = accelerator.Allocate1D(chunkCapacity); using MemoryBuffer1D resultBuffer = accelerator.Allocate1D(chunkCapacity); + using MemoryBuffer1D statusBuffer = accelerator.Allocate1D(chunkCapacity); int offset = 0; while (offset < length) { int batchSize = Math.Min(chunkCapacity, length - offset); Span divisorsSpan = divisors.AsSpan(0, batchSize); + Span destinationSpan = destination.AsSpan(offset, batchSize); + Span powSpan = powScratch.AsSpan(0, batchSize); + Span orderSpan = orderScratch.AsSpan(0, batchSize); + Span statusSpan = statusScratch.AsSpan(0, batchSize); + + int pending = 0; for (int i = 0; i < batchSize; i++) { - divisorsSpan[i] = start + (ulong)(offset + i); + ulong divisorValue = start + (ulong)(offset + i); + divisorsSpan[i] = divisorValue; + + if ((divisorValue & (divisorValue - 1UL)) == 0UL) + { + destinationSpan[i] = 1UL; + powSpan[i] = 0UL; + orderSpan[i] = 1UL; + statusSpan[i] = ByteOne; + } + else + { + ulong initialPow = 2UL; + if (initialPow >= divisorValue) + { + initialPow -= divisorValue; + } + + destinationSpan[i] = 0UL; + powSpan[i] = initialPow; + orderSpan[i] = 1UL; + statusSpan[i] = ByteZero; + pending++; + } } ArrayView1D divisorView = divisorBuffer.View.SubView(0, batchSize); + ArrayView1D powView = powBuffer.View.SubView(0, batchSize); + ArrayView1D orderView = orderBuffer.View.SubView(0, batchSize); ArrayView1D resultView = resultBuffer.View.SubView(0, batchSize); + ArrayView1D statusView = statusBuffer.View.SubView(0, batchSize); divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorsSpan), batchSize); - kernel(batchSize, divisorView, resultView); - accelerator.Synchronize(); + powView.CopyFromCPU(ref MemoryMarshal.GetReference(powSpan), batchSize); + orderView.CopyFromCPU(ref MemoryMarshal.GetReference(orderSpan), batchSize); + resultView.CopyFromCPU(ref MemoryMarshal.GetReference(destinationSpan), batchSize); + statusView.CopyFromCPU(ref MemoryMarshal.GetReference(statusSpan), batchSize); + + while (pending > 0) + { + kernel(batchSize, divisorView, powView, orderView, resultView, statusView); + accelerator.Synchronize(); + + statusView.CopyToCPU(ref MemoryMarshal.GetReference(statusSpan), batchSize); + + pending = 0; + for (int i = 0; i < batchSize; i++) + { + if (statusSpan[i] == ByteZero) + { + pending++; + } + } + } - Span destinationSpan = destination.AsSpan(offset, batchSize); resultView.CopyToCPU(ref MemoryMarshal.GetReference(destinationSpan), batchSize); + offset += batchSize; } } finally { ArrayPool.Shared.Return(divisors, clearArray: false); + ArrayPool.Shared.Return(powScratch, clearArray: false); + ArrayPool.Shared.Return(orderScratch, clearArray: false); + ArrayPool.Shared.Return(statusScratch, clearArray: false); } } finally @@ -373,14 +436,52 @@ public sealed class DivisorCycleCache } } - private static Action, ArrayView1D> LoadKernel(Accelerator accelerator) + private static Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> LoadKernel(Accelerator accelerator) { - return accelerator.LoadAutoGroupedStreamKernel, ArrayView1D>(GpuDivisorCycleKernel); + return accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(GpuAdvanceDivisorCyclesKernel); } - private static void GpuDivisorCycleKernel(Index1D index, ArrayView1D divisors, ArrayView1D cycles) + private static void GpuAdvanceDivisorCyclesKernel( + Index1D index, + ArrayView1D divisors, + ArrayView1D pow, + ArrayView1D order, + ArrayView1D cycles, + ArrayView1D status) { - cycles[index] = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[index]); + if (status[index] != ByteZero) + { + return; + } + + ulong divisor = divisors[index]; + ulong currentPow = pow[index]; + ulong currentOrder = order[index]; + int steps = GpuCycleStepsPerInvocation; + + do + { + currentPow += currentPow; + if (currentPow >= divisor) + { + currentPow -= divisor; + } + + currentOrder++; + + if (currentPow == 1UL) + { + cycles[index] = currentOrder; + status[index] = ByteOne; + pow[index] = currentPow; + order[index] = currentOrder; + return; + } + } + while (--steps != 0); + + pow[index] = currentPow; + order[index] = currentOrder; } private sealed class AcceleratorReferenceComparer : IEqualityComparer