diff --git a/lib/simulator_cuda.h b/lib/simulator_cuda.h index 5743bea8b..2a5e651e4 100644 --- a/lib/simulator_cuda.h +++ b/lib/simulator_cuda.h @@ -344,13 +344,13 @@ class SimulatorCUDA final { unsigned k = 5 + G; unsigned n = num_qubits > k ? num_qubits - k : 0; - unsigned size = unsigned{1} << n; + uint64_t size = uint64_t{1} << n; unsigned threads = 64U; - unsigned blocks = std::max(1U, size / 2); + uint64_t blocks = std::max(uint64_t{1}, size / 2); IndicesH d_i(d_ws); - ApplyGateH_Kernel<<>>( + ApplyGateH_Kernel<<>>( (fp_type*) d_ws, d_i.xss, d_i.ms, state.get()); } @@ -368,13 +368,13 @@ class SimulatorCUDA final { unsigned k = 5 + num_effective_qs; unsigned n = num_qubits > k ? num_qubits - k : 0; - unsigned size = unsigned{1} << n; + uint64_t size = uint64_t{1} << n; unsigned threads = 32; - unsigned blocks = size; + uint64_t blocks = size; IndicesL d_i(d_ws); - ApplyGateL_Kernel<<>>( + ApplyGateL_Kernel<<>>( (fp_type*) d_ws, d_i.xss, d_i.ms, d_i.qis, d_i.tis, 1 << num_effective_qs, state.get()); } @@ -401,13 +401,13 @@ class SimulatorCUDA final { unsigned k = 5 + G + cqs.size(); unsigned n = num_qubits > k ? num_qubits - k : 0; - unsigned size = unsigned{1} << n; + uint64_t size = uint64_t{1} << n; unsigned threads = 64U; - unsigned blocks = std::max(1U, size / 2); + uint64_t blocks = std::max(uint64_t{1}, size / 2); IndicesH d_i(d_ws); - ApplyControlledGateH_Kernel<<>>( + ApplyControlledGateH_Kernel<<>>( (fp_type*) d_ws, d_i.xss, d_i.ms, num_aqs + 1, cvalsh, state.get()); } @@ -426,13 +426,13 @@ class SimulatorCUDA final { unsigned k = 5 + G + cqs.size(); unsigned n = num_qubits > k ? num_qubits - k : 0; - unsigned size = unsigned{1} << n; + uint64_t size = uint64_t{1} << n; unsigned threads = 32; - unsigned blocks = size; + uint64_t blocks = size; IndicesL d_i(d_ws); - ApplyControlledGateLH_Kernel<<>>( + ApplyControlledGateLH_Kernel<<>>( (fp_type*) d_ws, d_i.xss, d_i.ms, d_i.qis, d_i.tis, d.num_aqs + 1, d.cvalsh, 1 << d.num_effective_qs, state.get()); } @@ -452,13 +452,13 @@ class SimulatorCUDA final { unsigned k = 5 + G + cqs.size(); unsigned n = num_qubits > k ? num_qubits - k : 0; - unsigned size = unsigned{1} << n; + uint64_t size = uint64_t{1} << n; unsigned threads = 32; - unsigned blocks = size; + uint64_t blocks = size; IndicesLC d_i(d_ws); - ApplyControlledGateL_Kernel<<>>( + ApplyControlledGateL_Kernel<<>>( (fp_type*) d_ws, d_i.xss, d_i.ms, d_i.qis, d_i.tis, d_i.cis, d.num_aqs + 1, d.cvalsh, 1 << d.num_effective_qs, 1 << (5 - d.remaining_low_cqs), state.get()); @@ -479,11 +479,11 @@ class SimulatorCUDA final { unsigned k = 5 + G; unsigned n = num_qubits > k ? num_qubits - k : 0; - unsigned size = unsigned{1} << n; + uint64_t size = uint64_t{1} << n; unsigned s = std::min(n >= 14 ? n - 14 : 0, 4U); unsigned threads = 64U; - unsigned blocks = std::max(1U, (size / 2) >> s); + uint64_t blocks = std::max(uint64_t{1}, (size / 2) >> s); unsigned num_iterations_per_block = 1 << s; constexpr unsigned m = 16; @@ -493,7 +493,7 @@ class SimulatorCUDA final { IndicesH d_i(d_ws); - ExpectationValueH_Kernel<<>>( + ExpectationValueH_Kernel<<>>( (fp_type*) d_ws, d_i.xss, d_i.ms, num_iterations_per_block, state.get(), Plus(), d_res1); @@ -517,11 +517,11 @@ class SimulatorCUDA final { unsigned k = 5 + num_effective_qs; unsigned n = num_qubits > k ? num_qubits - k : 0; - unsigned size = unsigned{1} << n; + uint64_t size = uint64_t{1} << n; unsigned s = std::min(n >= 13 ? n - 13 : 0, 5U); unsigned threads = 32; - unsigned blocks = size >> s; + uint64_t blocks = size >> s; unsigned num_iterations_per_block = 1 << s; constexpr unsigned m = 16; @@ -531,7 +531,7 @@ class SimulatorCUDA final { IndicesL d_i(d_ws); - ExpectationValueL_Kernel<<>>( + ExpectationValueL_Kernel<<>>( (fp_type*) d_ws, d_i.xss, d_i.ms, d_i.qis, d_i.tis, num_iterations_per_block, state.get(), Plus(), d_res1); @@ -542,7 +542,7 @@ class SimulatorCUDA final { template std::complex ExpectationValueReduceFinal( - unsigned blocks, double mul, + uint64_t blocks, double mul, const Complex* d_res1, Complex* d_res2) const { Complex res2[m]; @@ -550,10 +550,10 @@ class SimulatorCUDA final { ErrorCheck(cudaMemcpy(res2, d_res1, blocks * sizeof(Complex), cudaMemcpyDeviceToHost)); } else { - unsigned threads2 = std::min(1024U, blocks); - unsigned blocks2 = std::min(m, blocks / threads2); + unsigned threads2 = std::min(uint64_t{1024}, blocks); + uint64_t blocks2 = std::min(uint64_t{m}, blocks / threads2); - unsigned dblocks = std::max(1U, blocks / (blocks2 * threads2)); + unsigned dblocks = std::max(uint64_t{1}, blocks / (blocks2 * threads2)); unsigned bytes = threads2 * sizeof(Complex); Reduce2Kernel<<>>( @@ -568,7 +568,7 @@ class SimulatorCUDA final { double re = 0; double im = 0; - for (unsigned i = 0; i < blocks; ++i) { + for (uint64_t i = 0; i < blocks; ++i) { re += res2[i].re; im += res2[i].im; } diff --git a/lib/simulator_cuda_kernels.h b/lib/simulator_cuda_kernels.h index e21a9d62e..d39661c56 100644 --- a/lib/simulator_cuda_kernels.h +++ b/lib/simulator_cuda_kernels.h @@ -18,13 +18,13 @@ #ifdef __NVCC__ #include #include - - #include "util_cuda.h" #elif __HIP__ #include #include "cuda2hip.h" #endif +#include "util_cuda.h" + namespace qsim { template @@ -33,6 +33,8 @@ __global__ void ApplyGateH_Kernel( const idx_type* __restrict__ mss, fp_type* __restrict__ rstate) { // blockDim.x must be equal to 64. + uint64_t blockId = GetBlockId(); + static_assert(G < 7, "gates acting on more than 6 qubits are not supported."); constexpr unsigned gsize = 1 << G; @@ -61,7 +63,7 @@ __global__ void ApplyGateH_Kernel( __syncthreads(); - idx_type i = (64 * idx_type{blockIdx.x} + threadIdx.x) & 0xffffffffffe0; + idx_type i = (64 * idx_type{blockId} + threadIdx.x) & 0xffffffffffe0; idx_type ii = i & mss[0]; for (unsigned j = 1; j <= G; ++j) { i *= 2; @@ -115,6 +117,8 @@ __global__ void ApplyGateL_Kernel( fp_type* __restrict__ rstate) { // blockDim.x must be equal to 32. + uint64_t blockId = GetBlockId(); + static_assert(G < 7, "gates acting on more than 6 qubits are not supported."); constexpr unsigned gsize = 1 << G; @@ -137,7 +141,7 @@ __global__ void ApplyGateL_Kernel( } } - idx_type i = 32 * idx_type{blockIdx.x}; + idx_type i = 32 * idx_type{blockId}; idx_type ii = i & mss[0]; for (unsigned j = 1; j <= G; ++j) { i *= 2; @@ -204,6 +208,8 @@ __global__ void ApplyControlledGateH_Kernel( fp_type* __restrict__ rstate) { // blockDim.x must be equal to 64. + uint64_t blockId = GetBlockId(); + static_assert(G < 7, "gates acting on more than 6 qubits are not supported."); constexpr unsigned gsize = 1 << G; @@ -232,7 +238,7 @@ __global__ void ApplyControlledGateH_Kernel( __syncthreads(); - idx_type i = (64 * idx_type{blockIdx.x} + threadIdx.x) & 0xffffffffffe0; + idx_type i = (64 * idx_type{blockId} + threadIdx.x) & 0xffffffffffe0; idx_type ii = i & mss[0]; for (unsigned j = 1; j < num_mss; ++j) { i *= 2; @@ -288,6 +294,8 @@ __global__ void ApplyControlledGateLH_Kernel( unsigned esize, fp_type* __restrict__ rstate) { // blockDim.x must be equal to 32. + uint64_t blockId = GetBlockId(); + static_assert(G < 7, "gates acting on more than 6 qubits are not supported."); constexpr unsigned gsize = 1 << G; @@ -300,7 +308,7 @@ __global__ void ApplyControlledGateLH_Kernel( __shared__ fp_type rs0[32][gsize + 1], is0[32][gsize + 1]; __shared__ fp_type v[2 * gsize * rows]; - idx_type i = 32 * idx_type{blockIdx.x}; + idx_type i = 32 * idx_type{blockId}; idx_type ii = i & mss[0]; for (unsigned j = 1; j < num_mss; ++j) { i *= 2; @@ -381,6 +389,8 @@ __global__ void ApplyControlledGateL_Kernel( fp_type* __restrict__ rstate) { // blockDim.x must be equal to 32. + uint64_t blockId = GetBlockId(); + static_assert(G < 7, "gates acting on more than 6 qubits are not supported."); constexpr unsigned gsize = 1 << G; @@ -393,7 +403,7 @@ __global__ void ApplyControlledGateL_Kernel( __shared__ fp_type rs0[32][gsize + 1], is0[32][gsize + 1]; __shared__ fp_type v[2 * gsize * rows]; - idx_type i = 32 * idx_type{blockIdx.x}; + idx_type i = 32 * idx_type{blockId}; idx_type ii = i & mss[0]; for (unsigned j = 1; j < num_mss; ++j) { i *= 2; @@ -477,6 +487,8 @@ __global__ void ExpectationValueH_Kernel( const fp_type* __restrict__ rstate, Op op, cfp_type* __restrict__ result) { // blockDim.x must be equal to 64. + uint64_t blockId = GetBlockId(); + static_assert(G < 7, "gates acting on more than 6 qubits are not supported."); constexpr unsigned gsize = 1 << G; @@ -508,7 +520,7 @@ __global__ void ExpectationValueH_Kernel( double im = 0; for (unsigned iter = 0; iter < num_iterations_per_block; ++iter) { - idx_type b = num_iterations_per_block * idx_type{blockIdx.x} + iter; + idx_type b = num_iterations_per_block * idx_type{blockId} + iter; idx_type i = (64 * b + threadIdx.x) & 0xffffffffffe0; idx_type ii = i & mss[0]; @@ -573,8 +585,8 @@ __global__ void ExpectationValueH_Kernel( __syncthreads(); if (threadIdx.x == 0) { - result[blockIdx.x].re = partial2[0].re + partial2[1].re; - result[blockIdx.x].im = partial2[0].im + partial2[1].im; + result[blockId].re = partial2[0].re + partial2[1].re; + result[blockId].im = partial2[0].im + partial2[1].im; } } @@ -587,6 +599,8 @@ __global__ void ExpectationValueL_Kernel( const fp_type* __restrict__ rstate, Op op, cfp_type* __restrict__ result) { // blockDim.x must be equal to 32. + uint64_t blockId = GetBlockId(); + static_assert(G < 7, "gates acting on more than 6 qubits are not supported."); constexpr unsigned gsize = 1 << G; @@ -612,7 +626,7 @@ __global__ void ExpectationValueL_Kernel( double im = 0; for (idx_type iter = 0; iter < num_iterations_per_block; ++iter) { - idx_type i = 32 * (num_iterations_per_block * idx_type{blockIdx.x} + iter); + idx_type i = 32 * (num_iterations_per_block * idx_type{blockId} + iter); idx_type ii = i & mss[0]; for (unsigned j = 1; j <= G; ++j) { i *= 2; @@ -673,8 +687,8 @@ __global__ void ExpectationValueL_Kernel( auto val = WarpReduce(partial[threadIdx.x], op); if (threadIdx.x == 0) { - result[blockIdx.x].re = val.re; - result[blockIdx.x].im = val.im; + result[blockId].re = val.re; + result[blockId].im = val.im; } } diff --git a/lib/statespace_cuda.h b/lib/statespace_cuda.h index 660db074c..e0f3cfe03 100644 --- a/lib/statespace_cuda.h +++ b/lib/statespace_cuda.h @@ -49,7 +49,7 @@ class StateSpaceCUDA : struct Grid { unsigned threads; unsigned dblocks; - unsigned blocks; + dim3 blocks; }; public: @@ -86,10 +86,10 @@ class StateSpaceCUDA : uint64_t size = MinSize(state.num_qubits()) / 2; unsigned threads = std::min(size, uint64_t{param_.num_threads}); - unsigned blocks = size / threads; + uint64_t blocks = size / threads; unsigned bytes = 2 * threads * sizeof(fp_type); - InternalToNormalOrderKernel<<>>(state.get()); + InternalToNormalOrderKernel<<>>(state.get()); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); } @@ -98,10 +98,10 @@ class StateSpaceCUDA : uint64_t size = MinSize(state.num_qubits()) / 2; unsigned threads = std::min(size, uint64_t{param_.num_threads}); - unsigned blocks = size / threads; + uint64_t blocks = size / threads; unsigned bytes = 2 * threads * sizeof(fp_type); - NormalToInternalOrderKernel<<>>(state.get()); + NormalToInternalOrderKernel<<>>(state.get()); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); } @@ -117,11 +117,11 @@ class StateSpaceCUDA : uint64_t hsize = uint64_t{1} << state.num_qubits(); unsigned threads = std::min(size, uint64_t{param_.num_threads}); - unsigned blocks = size / threads; + uint64_t blocks = size / threads; fp_type v = double{1} / std::sqrt(hsize); - SetStateUniformKernel<<>>(v, hsize, state.get()); + SetStateUniformKernel<<>>(v, hsize, state.get()); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); } @@ -178,9 +178,9 @@ class StateSpaceCUDA : uint64_t size = MinSize(state.num_qubits()) / 2; unsigned threads = std::min(size, uint64_t{param_.num_threads}); - unsigned blocks = size / threads; + uint64_t blocks = size / threads; - BulkSetAmplKernel<<>>( + BulkSetAmplKernel<<>>( mask, bits, re, im, exclude, state.get()); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); @@ -195,9 +195,9 @@ class StateSpaceCUDA : uint64_t size = MinSize(src.num_qubits()); unsigned threads = std::min(size, uint64_t{param_.num_threads}); - unsigned blocks = size / threads; + uint64_t blocks = size / threads; - AddKernel<<>>(src.get(), dest.get()); + AddKernel<<>>(src.get(), dest.get()); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); @@ -209,9 +209,9 @@ class StateSpaceCUDA : uint64_t size = MinSize(state.num_qubits()); unsigned threads = std::min(size, uint64_t{param_.num_threads}); - unsigned blocks = size / threads; + uint64_t blocks = size / threads; - MultiplyKernel<<>>(a, state.get()); + MultiplyKernel<<>>(a, state.get()); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); } @@ -248,15 +248,17 @@ class StateSpaceCUDA : if (num_samples > 0) { Grid g1 = GetGrid1(MinSize(state.num_qubits()) / 2); unsigned bytes = g1.threads * sizeof(double); + uint64_t num_blocks1 = + MinSize(state.num_qubits()) / 2 / (g1.threads * g1.dblocks); - unsigned scratch_size = (g1.blocks + 1) * sizeof(double) + unsigned scratch_size = (num_blocks1 + 1) * sizeof(double) + num_samples * (sizeof(uint64_t) + sizeof(DistrRealType)); void* scratch = AllocScratch(scratch_size); double* d_res2 = (double*) scratch; double* d_res1 = d_res2 + 1; - uint64_t* d_bitstrings = (uint64_t*) (d_res1 + g1.blocks); + uint64_t* d_bitstrings = (uint64_t*) (d_res1 + num_blocks1); DistrRealType* d_rs = (DistrRealType *) (d_bitstrings + num_samples); auto op1 = RealProduct(); @@ -269,17 +271,17 @@ class StateSpaceCUDA : double norm; - if (g1.blocks == 1) { + if (num_blocks1 == 1) { ErrorCheck( cudaMemcpy(&norm, d_res1, sizeof(double), cudaMemcpyDeviceToHost)); } else { - Grid g2 = GetGrid2(g1.blocks); + Grid g2 = GetGrid2(num_blocks1); unsigned bytes = g2.threads * sizeof(double); auto op3 = Plus(); Reduce2Kernel<<>>( - g2.dblocks, g1.blocks, op3, op3, d_res1, d_res2); + g2.dblocks, num_blocks1, op3, op3, d_res1, d_res2); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); @@ -294,7 +296,7 @@ class StateSpaceCUDA : num_samples * sizeof(DistrRealType), cudaMemcpyHostToDevice)); - SampleKernel<<<1, g1.threads>>>(g1.blocks, g1.dblocks, num_samples, + SampleKernel<<<1, g1.threads>>>(num_blocks1, g1.dblocks, num_samples, d_rs, d_res1, state.get(), d_bitstrings); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); @@ -319,17 +321,19 @@ class StateSpaceCUDA : uint64_t size = MinSize(state.num_qubits()) / 2; unsigned threads = std::min(size, uint64_t{param_.num_threads}); - unsigned blocks = size / threads; + uint64_t blocks = size / threads; - CollapseKernel<<>>(mr.mask, mr.bits, renorm, state.get()); + CollapseKernel<<>>(mr.mask, mr.bits, renorm, state.get()); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); } std::vector PartialNorms(const State& state) const { Grid g = GetGrid1(MinSize(state.num_qubits()) / 2); + uint64_t num_blocks = + MinSize(state.num_qubits()) / 2 / (g.threads * g.dblocks); - unsigned scratch_size = g.blocks * sizeof(double); + unsigned scratch_size = num_blocks * sizeof(double); unsigned bytes = g.threads * sizeof(double); double* d_res = (double*) AllocScratch(scratch_size); @@ -342,7 +346,7 @@ class StateSpaceCUDA : ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); - std::vector norms(g.blocks); + std::vector norms(num_blocks); ErrorCheck( cudaMemcpy(norms.data(), d_res, scratch_size, cudaMemcpyDeviceToHost)); @@ -387,21 +391,18 @@ class StateSpaceCUDA : Grid GetGrid1(uint64_t size) const { Grid grid; - grid.threads = std::min(size, uint64_t{param_.num_threads}); grid.dblocks = std::min(size / grid.threads, uint64_t{param_.num_dblocks}); - grid.blocks = size / (grid.threads * grid.dblocks); - + uint64_t num_blocks = size / (grid.threads * grid.dblocks); + grid.blocks = CreateGrid(num_blocks); return grid; } Grid GetGrid2(unsigned size) const { Grid grid; - grid.threads = std::min(param_.num_threads, std::max(32U, size)); grid.dblocks = std::max(1U, size / grid.threads); - grid.blocks = 1; - + grid.blocks = dim3(1, 1, 1); return grid; } @@ -417,8 +418,9 @@ class StateSpaceCUDA : Grid g1 = GetGrid1(size); unsigned bytes = g1.threads * sizeof(FP1); + uint64_t num_blocks1 = size / (g1.threads * g1.dblocks); - FP2* d_res2 = (FP2*) AllocScratch((g1.blocks + 1) * sizeof(FP2)); + FP2* d_res2 = (FP2*) AllocScratch((num_blocks1 + 1) * sizeof(FP2)); FP2* d_res1 = d_res2 + 1; auto op1 = Op(); @@ -438,18 +440,18 @@ class StateSpaceCUDA : FP2 result; - if (g1.blocks == 1) { + if (num_blocks1 == 1) { ErrorCheck( cudaMemcpy(&result, d_res1, sizeof(FP2), cudaMemcpyDeviceToHost)); } else { - Grid g2 = GetGrid2(g1.blocks); + Grid g2 = GetGrid2(num_blocks1); unsigned bytes = g2.threads * sizeof(FP2); auto op2 = Plus(); auto op3 = Plus::type>(); Reduce2Kernel<<>>( - g2.dblocks, g1.blocks, op2, op3, d_res1, d_res2); + g2.dblocks, num_blocks1, op2, op3, d_res1, d_res2); ErrorCheck(cudaPeekAtLastError()); ErrorCheck(cudaDeviceSynchronize()); diff --git a/lib/statespace_cuda_kernels.h b/lib/statespace_cuda_kernels.h index 0bc4ba706..c762dfbb0 100644 --- a/lib/statespace_cuda_kernels.h +++ b/lib/statespace_cuda_kernels.h @@ -39,7 +39,7 @@ __device__ __forceinline__ FP1 BlockReduce1( unsigned warp = threadIdx.x / warp_size; unsigned lane = threadIdx.x % warp_size; - uint64_t k0 = 2 * n * blockIdx.x * blockDim.x + 2 * tid - lane; + uint64_t k0 = 2 * n * GetBlockId() * blockDim.x + 2 * tid - lane; uint64_t k1 = k0 + 2 * n * blockDim.x; FP1 r; @@ -88,7 +88,7 @@ __device__ __forceinline__ FP1 BlockReduce1Masked( unsigned warp = threadIdx.x / warp_size; unsigned lane = threadIdx.x % warp_size; - uint64_t k0 = 2 * n * blockIdx.x * blockDim.x + 2 * tid - lane; + uint64_t k0 = 2 * n * GetBlockId() * blockDim.x + 2 * tid - lane; uint64_t k1 = k0 + 2 * n * blockDim.x; FP1 r = 0; @@ -137,7 +137,7 @@ __device__ __forceinline__ FP1 BlockReduce2( FP1* partial1 = (FP1*) shared; unsigned tid = threadIdx.x; - uint64_t k0 = n * blockIdx.x * blockDim.x + tid; + uint64_t k0 = n * GetBlockId() * blockDim.x + tid; uint64_t k1 = k0 + n * blockDim.x; FP1 r = 0; @@ -185,7 +185,7 @@ __global__ void Reduce1Kernel(uint64_t n, Op1 op1, Op2 op2, Op3 op3, FP1 sum = detail::BlockReduce1(n, op1, op2, op3, s1, s2); if (threadIdx.x == 0) { - result[blockIdx.x] = sum; + result[GetBlockId()] = sum; } } @@ -198,7 +198,7 @@ __global__ void Reduce1MaskedKernel(uint64_t n, uint64_t mask, uint64_t bits, detail::BlockReduce1Masked(n, mask, bits, op1, op2, op3, s1, s2); if (threadIdx.x == 0) { - result[blockIdx.x] = sum; + result[GetBlockId()] = sum; } } @@ -209,7 +209,7 @@ __global__ void Reduce2Kernel( FP1 sum = detail::BlockReduce2(n, size, op2, op3, s); if (threadIdx.x == 0) { - result[blockIdx.x] = sum; + result[GetBlockId()] = sum; } } @@ -217,7 +217,7 @@ template __global__ void InternalToNormalOrderKernel(FP* state) { unsigned lane = threadIdx.x % warp_size; unsigned l = 2 * threadIdx.x - lane; - uint64_t k = 2 * uint64_t{blockIdx.x} * blockDim.x + l; + uint64_t k = 2 * GetBlockId() * blockDim.x + l; extern __shared__ float shared[]; FP* buf = (FP*) shared; @@ -235,7 +235,7 @@ template __global__ void NormalToInternalOrderKernel(FP* state) { unsigned lane = threadIdx.x % warp_size; unsigned l = 2 * threadIdx.x - lane; - uint64_t k = 2 * uint64_t{blockIdx.x} * blockDim.x + l; + uint64_t k = 2 * GetBlockId() * blockDim.x + l; extern __shared__ float shared[]; FP* buf = (FP*) shared; @@ -252,7 +252,7 @@ __global__ void NormalToInternalOrderKernel(FP* state) { template __global__ void SetStateUniformKernel(FP v, uint64_t size, FP* state) { unsigned lane = threadIdx.x % warp_size; - uint64_t k = 2 * (uint64_t{blockIdx.x} * blockDim.x + threadIdx.x) - lane; + uint64_t k = 2 * (GetBlockId() * blockDim.x + threadIdx.x) - lane; state[k] = lane < size ? v : 0; state[k + warp_size] = 0; @@ -260,19 +260,19 @@ __global__ void SetStateUniformKernel(FP v, uint64_t size, FP* state) { template __global__ void AddKernel(const FP* state1, FP* state2) { - uint64_t k = uint64_t{blockIdx.x} * blockDim.x + threadIdx.x; + uint64_t k = GetBlockId() * blockDim.x + threadIdx.x; state2[k] += state1[k]; } template __global__ void MultiplyKernel(FP a, FP* state) { - uint64_t k = uint64_t{blockIdx.x} * blockDim.x + threadIdx.x; + uint64_t k = GetBlockId() * blockDim.x + threadIdx.x; state[k] *= a; } template __global__ void CollapseKernel(uint64_t mask, uint64_t bits, FP r, FP* state) { - uint64_t k1 = uint64_t{blockIdx.x} * blockDim.x + threadIdx.x; + uint64_t k1 = GetBlockId() * blockDim.x + threadIdx.x; uint64_t k2 = 2 * k1 - threadIdx.x % warp_size; if ((k1 & mask) == bits) { @@ -287,7 +287,7 @@ __global__ void CollapseKernel(uint64_t mask, uint64_t bits, FP r, FP* state) { template __global__ void BulkSetAmplKernel( uint64_t mask, uint64_t bits, FP re, FP im, bool exclude, FP* state) { - uint64_t k1 = uint64_t{blockIdx.x} * blockDim.x + threadIdx.x; + uint64_t k1 = GetBlockId() * blockDim.x + threadIdx.x; uint64_t k2 = 2 * k1 - threadIdx.x % warp_size; bool set = ((k1 & mask) == bits) ^ exclude; @@ -299,7 +299,7 @@ __global__ void BulkSetAmplKernel( } template -__global__ void SampleKernel(unsigned num_blocks, +__global__ void SampleKernel(uint64_t num_blocks, uint64_t n, uint64_t num_samples, const FP1* rs, const FP2* ps, const FP3* state, uint64_t *bitstrings) { @@ -308,7 +308,7 @@ __global__ void SampleKernel(unsigned num_blocks, uint64_t m = 0; double csum = 0; - for (unsigned block_id = 0; block_id < num_blocks; ++block_id) { + for (uint64_t block_id = 0; block_id < num_blocks; ++block_id) { uint64_t km = n * blockDim.x; uint64_t k0 = block_id * km; diff --git a/lib/util_cuda.h b/lib/util_cuda.h index b34292753..004ddc11b 100644 --- a/lib/util_cuda.h +++ b/lib/util_cuda.h @@ -130,6 +130,25 @@ __device__ __forceinline__ FP1 WarpReduce(FP1 val, Op op) { return val; } +__device__ __forceinline__ uint64_t GetBlockId() { + return (uint64_t{blockIdx.z} * gridDim.y + uint64_t{blockIdx.y}) + * gridDim.x + blockIdx.x; +} + +inline dim3 CreateGrid(uint64_t blocks) { + if (blocks <= 65536) { + return dim3((uint32_t) blocks); + } + uint32_t x = 65536; + uint64_t rem = blocks / x; + if (rem <= 32768) { + return dim3(x, (uint32_t) rem); + } + uint32_t y = 32768; + uint32_t z = (uint32_t) (rem / y); + return dim3(x, y, z); +} + template __device__ __forceinline__ Complex WarpReduce(Complex val, Op op) { for (unsigned i = warp_size / 2; i > 0; i /= 2) { diff --git a/lib/vectorspace_cuda.h b/lib/vectorspace_cuda.h index 5cfd4e834..f1efdd51d 100644 --- a/lib/vectorspace_cuda.h +++ b/lib/vectorspace_cuda.h @@ -86,9 +86,8 @@ class VectorSpaceCUDA { static Vector Create(unsigned num_qubits) { fp_type* p; - auto size = sizeof(fp_type) * Impl::MinSize(num_qubits); + uint64_t size = sizeof(fp_type) * Impl::MinSize(num_qubits); auto rc = cudaMalloc(&p, size); - if (rc == cudaSuccess) { return Vector{Pointer{(fp_type*) p, &detail::free}, num_qubits}; } else { @@ -169,4 +168,4 @@ class VectorSpaceCUDA { } // namespace qsim -#endif // VECTORSPACE_CUDA_H_ +#endif // VECTORSPACE_CUDA_H_ \ No newline at end of file