diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs index db1c2fb..0cf9227 100644 --- a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs @@ -1,6 +1,3 @@ -using System; -using System.Collections.Generic; -using System.Threading; using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Jobs; using Open.Numeric.Primes; @@ -12,7 +9,6 @@ namespace EvenPerfectBitScanner.Benchmarks; [MemoryDiagnoser] public class PrimeTesterBenchmarks { - private static readonly HeuristicCombinedPrimeTester _heuristicCombinedTester = new(); private static readonly HeuristicPrimeTester _heuristicTester = new(); private ulong[] _candidates = []; @@ -34,7 +30,7 @@ public class PrimeTesterBenchmarks _ = Prime.Numbers.IsPrime(_candidates[0]); } - [Benchmark(Baseline = true)] + [Benchmark] public int HeuristicCpu() { HeuristicPrimeTester tester = _heuristicTester; @@ -104,18 +100,15 @@ public class PrimeTesterBenchmarks return primeCount; } - private static readonly PrimeTester _nonHeuristicTester = new(); - [Benchmark] public int NonHeuristicGpu() { - PrimeTester tester = _nonHeuristicTester; ulong[] values = _candidates; int primeCount = 0; for (int i = 0; i < values.Length; i++) { - if (tester.IsPrimeGpu(values[i])) + if (PrimeTester.IsPrimeGpu(values[i])) { primeCount++; } @@ -124,7 +117,7 @@ public class PrimeTesterBenchmarks return primeCount; } - [Benchmark] + [Benchmark(Baseline = true)] public int NonHeuristicCpu() { ulong[] values = _candidates; diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index bcd7937..48e8e27 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -12,7 +12,7 @@ namespace EvenPerfectBitScanner; internal static class Program { - private static ThreadLocal PrimeTesters = null!; + // private static ThreadLocal PrimeTesters = null!; private static ThreadLocal MersenneTesters = null!; private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP; private static MersenneNumberDivisorGpuTester? _divisorTester; @@ -176,7 +176,7 @@ internal static class Program blockSize = useByDivisor ? 1 : blockSize; CandidatesCalculator.Initialize(currentP); - PrimeTesters = new ThreadLocal(() => new PrimeTester(), trackAllValues: true); + // PrimeTesters = new ThreadLocal(() => new PrimeTester(), trackAllValues: true); // Note: --primes-device controls default device for library kernels; p primality remains CPU here. // Initialize per-thread p residue tracker (Identity model) at currentP if (!useDivisor && !_cliArguments.UseByDivisor) @@ -626,7 +626,7 @@ internal static class Program { if (!(_runPrimesOnCpu ? PrimeTester.IsPrime(p) - : PrimeTesters.Value!.IsPrimeGpu(p))) + : PrimeTester.IsPrimeGpu(p))) { _lastCompositeP = true; return false; diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs index 9ba37b0..8f6ed46 100644 --- a/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs @@ -16,7 +16,7 @@ public class PrimeTesterGpuTests ulong[] primes = [31UL, 61UL, 89UL, 107UL, 127UL, 521UL]; foreach (ulong prime in primes) { - tester.IsPrimeGpu(prime).Should().BeTrue(); + PrimeTester.IsPrimeGpu(prime).Should().BeTrue(); } } finally @@ -34,7 +34,7 @@ public class PrimeTesterGpuTests ulong[] composites = [33UL, 39UL, 51UL, 77UL, 91UL]; foreach (ulong composite in composites) { - tester.IsPrimeGpu(composite).Should().BeFalse(); + PrimeTester.IsPrimeGpu(composite).Should().BeFalse(); } } finally diff --git a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs index d1efe98..72dbcfb 100644 --- a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs +++ b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs @@ -39,7 +39,7 @@ public class PrimeTesterIsPrimeTests { var tester = new PrimeTester(); - tester.IsPrimeGpu(n).Should().Be(expected); + PrimeTester.IsPrimeGpu(n).Should().Be(expected); } [Fact] @@ -64,7 +64,7 @@ public class PrimeTesterIsPrimeTests var testers = new ThreadLocal(() => new PrimeTester(), trackAllValues: true); Parallel.For(0, primes.Length, i => { - results[i] = testers.Value!.IsPrimeGpu(primes[i]); + results[i] = PrimeTester.IsPrimeGpu(primes[i]); }); results.Should().AllBeEquivalentTo(true); @@ -73,7 +73,7 @@ public class PrimeTesterIsPrimeTests { foreach (ulong prime in primes) { - testers.Value!.IsPrimeGpu(prime).Should().BeTrue(); + PrimeTester.IsPrimeGpu(prime).Should().BeTrue(); } }); } diff --git a/PerfectNumbers.Core/Gpu/GpuContextPool.cs b/PerfectNumbers.Core/Gpu/GpuContextPool.cs index b6218c4..00cb47f 100644 --- a/PerfectNumbers.Core/Gpu/GpuContextPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuContextPool.cs @@ -1,4 +1,3 @@ -using System; using System.Collections.Concurrent; using ILGPU; using ILGPU.Runtime; @@ -15,14 +14,11 @@ public static class GpuContextPool public readonly Context Context; public readonly Accelerator Accelerator; - private static readonly Context _sharedContext = Context.CreateDefault(); - private static readonly Device _sharedDevice = _sharedContext.GetPreferredDevice(false); - private static readonly Accelerator _sharedAccelerator = _sharedDevice.CreateAccelerator(_sharedContext); - public PooledContext() { - Context = _sharedContext; - Accelerator = _sharedAccelerator; + Context = SharedGpuContext.Context; + Accelerator = SharedGpuContext.Accelerator; + GpuStaticTableInitializer.EnsureStaticTables(Accelerator); // NOTE: Avoid loading/compiling any kernel here to prevent implicit // CL stream/queue creation during accelerator construction. // Some OpenCL drivers are fragile when a queue is created immediately @@ -48,7 +44,7 @@ public static class GpuContextPool private static int WarmedGpuContextCount; // Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu. - public static GpuContextLease Rent() + public static GpuContextLease Rent() { if (GpuPool.TryDequeue(out var gpu)) { diff --git a/PerfectNumbers.Core/Gpu/GpuKernelLease.cs b/PerfectNumbers.Core/Gpu/GpuKernelLease.cs index a79af7a..9b84a9d 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelLease.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelLease.cs @@ -2,6 +2,7 @@ using System; using System.Collections.Concurrent; using ILGPU; using ILGPU.Runtime; +using PerfectNumbers.Core; using static PerfectNumbers.Core.Gpu.GpuContextPool; namespace PerfectNumbers.Core.Gpu; @@ -121,6 +122,20 @@ public sealed class GpuKernelLease } } + public Action, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D> SpecialMaxKernel + { + get + { + var accel = Accelerator; + return KernelContainer.InitOnce(ref _kernels!.SpecialMax, () => + { + var loaded = accel.LoadAutoGroupedStreamKernel, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>(PrimeOrderGpuHeuristics.EvaluateSpecialMaxCandidatesKernel); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>>(); + }); + } + } + public Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernel { get diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index 0257fc9..9e65ef3 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -1,3 +1,4 @@ +using PerfectNumbers.Core; using System.Collections.Concurrent; using ILGPU; using ILGPU.Runtime; @@ -37,6 +38,7 @@ public sealed class KernelContainer public Action, ArrayView1D>? Pow2ModOrder; public Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>? SmallPrimeFactor; + public Action, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>? SpecialMax; // Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length. public MemoryBuffer1D? SmallCycles; @@ -46,6 +48,9 @@ public sealed class KernelContainer public MemoryBuffer1D? SmallPrimeFactorExponentSlots; public MemoryBuffer1D? SmallPrimeFactorCountSlot; public MemoryBuffer1D? SmallPrimeFactorRemainingSlot; + public MemoryBuffer1D? SpecialMaxFactors; + public MemoryBuffer1D? SpecialMaxCandidates; + public MemoryBuffer1D? SpecialMaxResult; public MemoryBuffer1D? SmallPrimesLastOne; public MemoryBuffer1D? SmallPrimesPow2LastOne; public MemoryBuffer1D? SmallPrimesLastSeven; @@ -129,6 +134,29 @@ public readonly struct SmallPrimeFactorScratch } } +public readonly struct SpecialMaxScratch +{ + public readonly MemoryBuffer1D FactorValues; + public readonly MemoryBuffer1D CandidateValues; + public readonly MemoryBuffer1D ResultSlot; + + public SpecialMaxScratch( + MemoryBuffer1D factorValues, + MemoryBuffer1D candidateValues, + MemoryBuffer1D resultSlot) + { + FactorValues = factorValues; + CandidateValues = candidateValues; + ResultSlot = resultSlot; + } + + public ArrayView1D FactorsView => FactorValues.View; + + public ArrayView1D CandidatesView => CandidateValues.View; + + public ArrayView1D ResultView => ResultSlot.View; +} + public class GpuKernelPool { private static readonly ConcurrentDictionary KernelCache = new(); // TODO: Replace this concurrent map with a simple accelerator-indexed lookup once kernel launchers are prewarmed during startup so we can drop the thread-safe wrapper entirely. @@ -138,6 +166,13 @@ public class GpuKernelPool return KernelCache.GetOrAdd(accelerator, _ => new KernelContainer()); } + internal static void PreloadStaticTables(Accelerator accelerator) + { + EnsureSmallCyclesOnDevice(accelerator); + _ = EnsureSmallPrimesOnDevice(accelerator); + _ = EnsureSmallPrimeFactorTables(accelerator); + } + // Ensures the small cycles table is uploaded to the device for the given accelerator. // Returns the ArrayView to pass into kernels (when kernels are extended to accept it). public static ArrayView1D EnsureSmallCyclesOnDevice(Accelerator accelerator) @@ -277,6 +312,43 @@ public class GpuKernelPool } } + public static SpecialMaxScratch EnsureSpecialMaxScratch(Accelerator accelerator, int factorCapacity) + { + if (factorCapacity <= 0) + { + factorCapacity = 1; + } + + var kernels = GetKernels(accelerator); + lock (kernels) + { + MemoryBuffer1D? factorValues = kernels.SpecialMaxFactors; + if (factorValues is null || factorValues.Length < factorCapacity) + { + factorValues?.Dispose(); + factorValues = accelerator.Allocate1D(factorCapacity); + kernels.SpecialMaxFactors = factorValues; + } + + MemoryBuffer1D? candidateValues = kernels.SpecialMaxCandidates; + if (candidateValues is null || candidateValues.Length < factorCapacity) + { + candidateValues?.Dispose(); + candidateValues = accelerator.Allocate1D(factorCapacity); + kernels.SpecialMaxCandidates = candidateValues; + } + + MemoryBuffer1D? resultSlot = kernels.SpecialMaxResult; + if (resultSlot is null) + { + resultSlot = accelerator.Allocate1D(1); + kernels.SpecialMaxResult = resultSlot; + } + + return new SpecialMaxScratch(factorValues!, candidateValues!, resultSlot!); + } + } + public static GpuKernelLease GetKernel(bool useGpuOrder) { var limiter = GpuPrimeWorkLimiter.Acquire(); diff --git a/PerfectNumbers.Core/Gpu/GpuStaticTableInitializer.cs b/PerfectNumbers.Core/Gpu/GpuStaticTableInitializer.cs new file mode 100644 index 0000000..247ee09 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/GpuStaticTableInitializer.cs @@ -0,0 +1,44 @@ +using System.Collections.Concurrent; +using System.Collections.Generic; +using System.Runtime.CompilerServices; +using ILGPU; +using ILGPU.Runtime; +using PerfectNumbers.Core; + +namespace PerfectNumbers.Core.Gpu; + +internal static class GpuStaticTableInitializer +{ + private static readonly ConcurrentDictionary Initialized = new(AcceleratorReferenceComparer.Instance); + + internal static void EnsureStaticTables(Accelerator accelerator) + { + if (accelerator is null) + { + return; + } + + Initialized.GetOrAdd(accelerator, static acc => + { + PrimeTester.PrimeTesterGpuContextPool.EnsureStaticTables(acc); + GpuKernelPool.PreloadStaticTables(acc); + PrimeOrderGpuHeuristics.PreloadStaticTables(acc); + return 0; + }); + } + + private sealed class AcceleratorReferenceComparer : IEqualityComparer + { + internal static AcceleratorReferenceComparer Instance { get; } = new(); + + public bool Equals(Accelerator? x, Accelerator? y) + { + return ReferenceEquals(x, y); + } + + public int GetHashCode(Accelerator obj) + { + return RuntimeHelpers.GetHashCode(obj); + } + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs new file mode 100644 index 0000000..858e2a9 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs @@ -0,0 +1,80 @@ +using ILGPU; +using ILGPU.Runtime; +using PerfectNumbers.Core; + +namespace PerfectNumbers.Core.Gpu; + +internal static partial class PrimeOrderGpuHeuristics +{ + internal static void EvaluateSpecialMaxCandidatesKernel( + Index1D index, + ulong phi, + ArrayView1D factors, + int factorCount, + MontgomeryDivisorData divisor, + ArrayView1D candidates, + ArrayView1D resultOut) + { + if (index != 0) + { + return; + } + + int actual = 0; + for (int i = 0; i < factorCount; i++) + { + ulong factor = factors[i]; + if (factor <= 1UL) + { + continue; + } + + ulong reduced = phi / factor; + if (reduced == 0UL) + { + continue; + } + + candidates[actual] = reduced; + actual++; + } + + if (actual == 0) + { + resultOut[0] = 1; + return; + } + + for (int i = 1; i < actual; i++) + { + ulong current = candidates[i]; + int insert = i - 1; + while (insert >= 0 && candidates[insert] > current) + { + candidates[insert + 1] = candidates[insert]; + insert--; + } + + candidates[insert + 1] = current; + } + + var divisorPartial = new GpuDivisorPartialData(divisor.Modulus); + var stepper = new ExponentRemainderStepperGpu(divisorPartial); + if (stepper.InitializeIsUnityGpu(candidates[0])) + { + resultOut[0] = 0; + return; + } + + for (int i = 1; i < actual; i++) + { + if (stepper.ComputeNextIsUnityGpu(candidates[i])) + { + resultOut[0] = 0; + return; + } + } + + resultOut[0] = 1; + } +} diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index bfb9851..8664d8b 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -18,6 +18,11 @@ internal enum GpuPow2ModStatus internal static partial class PrimeOrderGpuHeuristics { + internal static void PreloadStaticTables(Accelerator accelerator) + { + _ = GetSmallPrimeDeviceCache(accelerator); + } + private static readonly ConcurrentDictionary OverflowedPrimes = new(); private static readonly ConcurrentDictionary OverflowedPrimesWide = new(); private static readonly ConcurrentDictionary, MontgomeryDivisorData, ArrayView1D>> Pow2ModKernelCache = new(); diff --git a/PerfectNumbers.Core/Gpu/SharedGpuContext.cs b/PerfectNumbers.Core/Gpu/SharedGpuContext.cs new file mode 100644 index 0000000..ec5c4cb --- /dev/null +++ b/PerfectNumbers.Core/Gpu/SharedGpuContext.cs @@ -0,0 +1,19 @@ +using System; +using System.Threading; +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal static class SharedGpuContext +{ + private static readonly Lazy s_context = new(() => ILGPU.Context.CreateDefault(), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy s_device = new(() => s_context.Value.GetPreferredDevice(false), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy s_accelerator = new(() => s_device.Value.CreateAccelerator(s_context.Value), LazyThreadSafetyMode.ExecutionAndPublication); + + internal static Context Context => s_context.Value; + + internal static Device Device => s_device.Value; + + internal static Accelerator Accelerator => s_accelerator.Value; +} diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 149bc77..9e7b026 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -412,7 +412,7 @@ public class MersenneDivisorCycles // Console.WriteLine($"MersenneDivisorCycles.TryFactorIntoCountsInternal hits {Volatile.Read(ref _tryFactorIntoCountsInternalHits)}"); // HeuristicPrimeTester primeTester = _primeTester ??= new(); - if (PrimeTester.IsPrime(remaining)) + if (PrimeTester.IsPrimeGpu(remaining)) // if (PrimeTester.IsPrimeCpu(remaining, CancellationToken.None)) { AddFactor(counts, remaining); @@ -943,7 +943,7 @@ public class MersenneDivisorCycles // HeuristicPrimeTester primeTester = _primeTester ??= new(); // if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrimeCpu(divisor, CancellationToken.None)) // if (!skipPrimeOrderHeuristic && primeTester.IsPrimeGpu(divisor)) - if (!skipPrimeOrderHeuristic && PrimeTester.IsPrime(divisor)) + if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeGpu(divisor)) { ulong computedOrder = PrimeOrderCalculator.Calculate( divisor, diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 702039e..49eb888 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -1,10 +1,8 @@ -using System; using System.Buffers; using System.Diagnostics; using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; -using ILGPU; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; @@ -108,7 +106,6 @@ internal static partial class PrimeOrderCalculator return candidateOrder; } - private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) { int length = factors.Count; @@ -121,16 +118,14 @@ internal static partial class PrimeOrderCalculator ReadOnlySpan factorSpan = new(factors.Factors, 0, length); - if (length <= 32) + Span stackBuffer = stackalloc ulong[length]; + + if (length <= 7) { - Span stackBuffer = stackalloc ulong[length]; return EvaluateSpecialMaxCandidates(stackBuffer, factorSpan, phi, prime, divisorData); } - ulong[] rented = ThreadStaticPools.UlongPool.Rent(length); - bool result = EvaluateSpecialMaxCandidates(rented.AsSpan(0, length), factorSpan, phi, prime, divisorData); - ThreadStaticPools.UlongPool.Return(rented, clearArray: false); - return result; + return EvaluateSpecialMaxCandidatesGpu(stackBuffer, factorSpan, phi, prime, divisorData); } private static bool EvaluateSpecialMaxCandidates(Span buffer, ReadOnlySpan factors, ulong phi, ulong prime, in MontgomeryDivisorData divisorData) @@ -1194,7 +1189,7 @@ internal static partial class PrimeOrderCalculator // Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}"); // bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); - bool isPrime = PrimeTester.IsPrime(composite); + bool isPrime = PrimeTester.IsPrimeGpu(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); entry = entry.WithPrimality(isPrime); @@ -1231,7 +1226,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorCofactorHits, 1UL); // Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}"); - cofactorIsPrime = PrimeTester.IsPrime(cofactor); + cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor); // cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor); // cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor); } @@ -1701,7 +1696,7 @@ internal static partial class PrimeOrderCalculator // HeuristicPrimeTester tester = _tester ??= new(); // bool isPrime = primeTester.HeuristicIsPrimeGpu(value); - bool isPrime = PrimeTester.IsPrime(value); + bool isPrime = PrimeTester.IsPrimeGpu(value); // bool isPrime = PrimeTester.IsPrimeGpu(value); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value); if (!knownComposite && isPrime) @@ -1717,7 +1712,7 @@ internal static partial class PrimeOrderCalculator // Console.WriteLine($"Factor completely after PollardRho hits {Volatile.Read(ref _factorCompletelyPollardRhoHits)}"); // isPrime = primeTester.HeuristicIsPrimeGpu(factor); - isPrime = PrimeTester.IsPrime(factor); + isPrime = PrimeTester.IsPrimeGpu(factor); // isPrime = PrimeTester.IsPrimeGpu(factor); // isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor); if (isPrime) @@ -1754,7 +1749,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorCofactorHits, 1UL); // Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}"); - isPrime = PrimeTester.IsPrime(quotient); + isPrime = PrimeTester.IsPrimeGpu(quotient); // isPrime = PrimeTester.IsPrimeGpu(quotient); // isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient); } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index 41ef94a..6dfd404 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -9,63 +9,112 @@ namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { - private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2; + private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2; - private const int GpuSmallPrimeFactorSlots = 64; + private const int GpuSmallPrimeFactorSlots = 64; - private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out int factorCount, out ulong remaining) - { - var primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); - var exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); - Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - remaining = value; + private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out int factorCount, out ulong remaining) + { + var primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); + var exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); + Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + remaining = value; - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); - Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; - SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator); - SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); + SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator); + SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); - var kernel = lease.SmallPrimeFactorKernel; - kernel( - stream, - 1, - value, - limit, - tables.PrimesView, - tables.SquaresView, - tables.Count, - scratch.PrimeSlots.View, - scratch.ExponentSlots.View, - scratch.CountSlot.View, - scratch.RemainingSlot.View); + var kernel = lease.SmallPrimeFactorKernel; + kernel( + stream, + 1, + value, + limit, + tables.PrimesView, + tables.SquaresView, + tables.Count, + scratch.PrimeSlots.View, + scratch.ExponentSlots.View, + scratch.CountSlot.View, + scratch.RemainingSlot.View); - stream.Synchronize(); + stream.Synchronize(); - factorCount = 0; - scratch.CountSlot.View.CopyToCPU(ref factorCount, 1); - factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots); + factorCount = 0; + scratch.CountSlot.View.CopyToCPU(ref factorCount, 1); + factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots); - if (factorCount > 0) - { - scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount); - scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount); - } + if (factorCount > 0) + { + scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount); + scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount); + } - scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1); + scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1); - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeBuffer[i]; - int exponent = exponentBuffer[i]; - counts.Add(primeValue, exponent); - } + for (int i = 0; i < factorCount; i++) + { + ulong primeValue = primeBuffer[i]; + int exponent = exponentBuffer[i]; + counts.Add(primeValue, exponent); + } - lease.Dispose(); - ThreadStaticPools.UlongPool.Return(primeBufferArray); - ThreadStaticPools.IntPool.Return(exponentBufferArray); - return true; - } + lease.Dispose(); + ThreadStaticPools.UlongPool.Return(primeBufferArray); + ThreadStaticPools.IntPool.Return(exponentBufferArray); + return true; + } + + private static bool EvaluateSpecialMaxCandidatesGpu( + Span buffer, + ReadOnlySpan factors, + ulong phi, + ulong prime, + in MontgomeryDivisorData divisorData) + { + if (factors.Length == 0) + { + return true; + } + + int factorCount = factors.Length; + Span factorSpan = buffer[..factorCount]; + for (int i = 0; i < factorCount; i++) + { + factorSpan[i] = factors[i].Value; + } + + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + + SpecialMaxScratch scratch = GpuKernelPool.EnsureSpecialMaxScratch(accelerator, factorCount); + + scratch.FactorsView.SubView(0, factorCount).CopyFromCPU(ref MemoryMarshal.GetReference(factorSpan), factorCount); + + var kernel = lease.SpecialMaxKernel; + kernel( + stream, + 1, + phi, + scratch.FactorsView, + factorCount, + divisorData, + scratch.CandidatesView, + scratch.ResultView); + + stream.Synchronize(); + + byte result = 0; + scratch.ResultView.CopyToCPU(ref result, 1); + + lease.Dispose(); + + return result != 0; + } } + diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 75f2b0d..4682952 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -93,50 +93,43 @@ public sealed class PrimeTester } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool IsPrimeGpu(ulong n) + public static bool IsPrimeGpu(ulong n) { - // bool belowGpuRange = n < 31UL; + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTesterGpuContextPool.Rent(1); + var accelerator = gpu.Accelerator; + var kernel = gpu.Kernel; - // if (belowGpuRange) - // { - return IsPrime(n); - // } + ulong value = n; + byte flag = 0; - // var limiter = GpuPrimeWorkLimiter.Acquire(); - // var gpu = PrimeTesterGpuContextPool.Rent(1); - // var accelerator = gpu.Accelerator; - // var kernel = gpu.Kernel; + var input = gpu.Input; + var output = gpu.Output; + var inputView = input.View; + var outputView = output.View; - // ulong value = n; - // byte flag = 0; + inputView.CopyFromCPU(ref value, 1); + kernel( + 1, + inputView, + gpu.DevicePrimesDefault.View, + gpu.DevicePrimesLastOne.View, + gpu.DevicePrimesLastSeven.View, + gpu.DevicePrimesLastThree.View, + gpu.DevicePrimesLastNine.View, + gpu.DevicePrimesPow2Default.View, + gpu.DevicePrimesPow2LastOne.View, + gpu.DevicePrimesPow2LastSeven.View, + gpu.DevicePrimesPow2LastThree.View, + gpu.DevicePrimesPow2LastNine.View, + outputView); + accelerator.Synchronize(); + outputView.CopyToCPU(ref flag, 1); - // var input = gpu.Input; - // var output = gpu.Output; - // var inputView = input.View; - // var outputView = output.View; + gpu.Dispose(); + limiter.Dispose(); - // inputView.CopyFromCPU(ref value, 1); - // kernel( - // 1, - // inputView, - // gpu.DevicePrimesDefault.View, - // gpu.DevicePrimesLastOne.View, - // gpu.DevicePrimesLastSeven.View, - // gpu.DevicePrimesLastThree.View, - // gpu.DevicePrimesLastNine.View, - // gpu.DevicePrimesPow2Default.View, - // gpu.DevicePrimesPow2LastOne.View, - // gpu.DevicePrimesPow2LastSeven.View, - // gpu.DevicePrimesPow2LastThree.View, - // gpu.DevicePrimesPow2LastNine.View, - // outputView); - // accelerator.Synchronize(); - // outputView.CopyToCPU(ref flag, 1); - - // gpu.Dispose(); - // limiter.Dispose(); - - // return flag != 0; + return flag != 0; } public static int GpuBatchSize { get; set; } = 262_144; @@ -221,23 +214,16 @@ public sealed class PrimeTester private static readonly object WarmUpLock = new(); private static int WarmedLeaseCount; private static readonly object SharedTablesLock = new(); - private static readonly Dictionary SharedTables = new(); - - private static string CreateAcceleratorKey(Accelerator accelerator) - { - return accelerator.AcceleratorType + ":" + accelerator.Name; - } + private static readonly Dictionary SharedTables = new(AcceleratorReferenceComparer.Instance); private static SharedHeuristicGpuTables RentSharedTables(Accelerator accelerator) { - string key = CreateAcceleratorKey(accelerator); - lock (SharedTablesLock) { - if (!SharedTables.TryGetValue(key, out var tables)) + if (!SharedTables.TryGetValue(accelerator, out var tables)) { - tables = new SharedHeuristicGpuTables(accelerator, key); - SharedTables.Add(key, tables); + tables = new SharedHeuristicGpuTables(accelerator); + SharedTables.Add(accelerator, tables); } tables.AddReference(); @@ -258,6 +244,19 @@ public sealed class PrimeTester } } + internal static void EnsureStaticTables(Accelerator accelerator) + { + if (accelerator is null) + { + return; + } + + var tables = RentSharedTables(accelerator); + var heuristicTables = tables.CreateHeuristicDivisorTables(); + HeuristicGpuDivisorTables.InitializeShared(in heuristicTables); + ReleaseSharedTables(tables); + } + internal static void WarmUp(int target, int minBufferCapacity) { if (target <= 0) @@ -272,6 +271,9 @@ public sealed class PrimeTester return; } + Accelerator accelerator = SharedGpuContext.Accelerator; + GpuStaticTableInitializer.EnsureStaticTables(accelerator); + int toCreate = target - WarmedLeaseCount; for (int i = 0; i < toCreate; i++) { @@ -319,6 +321,8 @@ public sealed class PrimeTester } SharedTables.Clear(); + var emptyTables = default(HeuristicGpuDivisorTables); + HeuristicGpuDivisorTables.InitializeShared(in emptyTables); } WarmedLeaseCount = 0; @@ -354,6 +358,20 @@ public sealed class PrimeTester WarmedLeaseCount = Pool.Count; } + + lock (SharedTablesLock) + { + if (SharedTables.Remove(accelerator, out var tables)) + { + tables.Dispose(); + } + + if (SharedTables.Count == 0) + { + var emptyTables = default(HeuristicGpuDivisorTables); + HeuristicGpuDivisorTables.InitializeShared(in emptyTables); + } + } } internal sealed class PrimeTesterGpuContextLease : IDisposable @@ -401,36 +419,17 @@ public sealed class PrimeTester public HeuristicGpuDivisorTables HeuristicGpuTables => _heuristicDivisorTables; - private static readonly Context _sharedContext = Context.CreateDefault(); - private static readonly Device _sharedDevice = _sharedContext.GetPreferredDevice(false); - private static readonly Accelerator _sharedAccelerator = _sharedDevice.CreateAccelerator(_sharedContext); - internal PrimeTesterGpuContextLease(int minBufferCapacity) { - AcceleratorContext = _sharedContext;// Context.CreateDefault(); - Accelerator = _sharedAccelerator; //AcceleratorContext.GetPreferredDevice(false).CreateAccelerator(AcceleratorContext); + AcceleratorContext = SharedGpuContext.Context; + Accelerator = SharedGpuContext.Accelerator; + GpuStaticTableInitializer.EnsureStaticTables(Accelerator); Kernel = Accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); HeuristicTrialDivisionKernel = Accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, HeuristicGpuDivisorTableKind, HeuristicGpuDivisorTables>(PrimeTesterKernels.HeuristicTrialDivisionKernel); _sharedTables = RentSharedTables(Accelerator); - _heuristicDivisorTables = new HeuristicGpuDivisorTables( - _sharedTables.HeuristicCombinedDivisorsEnding1.View, - _sharedTables.HeuristicCombinedDivisorsEnding3.View, - _sharedTables.HeuristicCombinedDivisorsEnding7.View, - _sharedTables.HeuristicCombinedDivisorsEnding9.View, - _sharedTables.HeuristicCombinedDivisorSquaresEnding1.View, - _sharedTables.HeuristicCombinedDivisorSquaresEnding3.View, - _sharedTables.HeuristicCombinedDivisorSquaresEnding7.View, - _sharedTables.HeuristicCombinedDivisorSquaresEnding9.View, - _sharedTables.HeuristicGroupADivisors.View, - _sharedTables.HeuristicGroupADivisorSquares.View, - _sharedTables.HeuristicGroupBDivisorsEnding1.View, - _sharedTables.HeuristicGroupBDivisorSquaresEnding1.View, - _sharedTables.HeuristicGroupBDivisorsEnding7.View, - _sharedTables.HeuristicGroupBDivisorSquaresEnding7.View, - _sharedTables.HeuristicGroupBDivisorsEnding9.View, - _sharedTables.HeuristicGroupBDivisorSquaresEnding9.View); + _heuristicDivisorTables = _sharedTables.CreateHeuristicDivisorTables(); HeuristicGpuDivisorTables.InitializeShared(in _heuristicDivisorTables); EnsureCapacity(minBufferCapacity); @@ -481,13 +480,27 @@ public sealed class PrimeTester } } + private sealed class AcceleratorReferenceComparer : IEqualityComparer + { + internal static AcceleratorReferenceComparer Instance { get; } = new(); + + public bool Equals(Accelerator? x, Accelerator? y) + { + return ReferenceEquals(x, y); + } + + public int GetHashCode(Accelerator obj) + { + return RuntimeHelpers.GetHashCode(obj); + } + } + private sealed class SharedHeuristicGpuTables { private int _referenceCount; - internal SharedHeuristicGpuTables(Accelerator accelerator, string key) + internal SharedHeuristicGpuTables(Accelerator accelerator) { - Key = key; var primesDefault = DivisorGenerator.SmallPrimes; DevicePrimesDefault = accelerator.Allocate1D(primesDefault.Length); @@ -568,7 +581,6 @@ public sealed class PrimeTester HeuristicCombinedDivisorSquaresEnding9 = CopySpanToDevice(accelerator, combinedSquaresEnding9); } - internal string Key { get; } internal MemoryBuffer1D DevicePrimesDefault { get; } internal MemoryBuffer1D DevicePrimesLastOne { get; } internal MemoryBuffer1D DevicePrimesLastSeven { get; } @@ -596,6 +608,27 @@ public sealed class PrimeTester internal MemoryBuffer1D HeuristicCombinedDivisorsEnding9 { get; } internal MemoryBuffer1D HeuristicCombinedDivisorSquaresEnding9 { get; } + internal HeuristicGpuDivisorTables CreateHeuristicDivisorTables() + { + return new HeuristicGpuDivisorTables( + HeuristicCombinedDivisorsEnding1.View, + HeuristicCombinedDivisorsEnding3.View, + HeuristicCombinedDivisorsEnding7.View, + HeuristicCombinedDivisorsEnding9.View, + HeuristicCombinedDivisorSquaresEnding1.View, + HeuristicCombinedDivisorSquaresEnding3.View, + HeuristicCombinedDivisorSquaresEnding7.View, + HeuristicCombinedDivisorSquaresEnding9.View, + HeuristicGroupADivisors.View, + HeuristicGroupADivisorSquares.View, + HeuristicGroupBDivisorsEnding1.View, + HeuristicGroupBDivisorSquaresEnding1.View, + HeuristicGroupBDivisorsEnding7.View, + HeuristicGroupBDivisorSquaresEnding7.View, + HeuristicGroupBDivisorsEnding9.View, + HeuristicGroupBDivisorSquaresEnding9.View); + } + internal void AddReference() { _referenceCount++;