cuda-samples/Samples/simpleVulkanMMAP/MonteCarloPi.cu

/* Copyright (c) 2020, NVIDIA CORPORATION. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *  * Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *  * Neither the name of NVIDIA CORPORATION nor the names of its
 *    contributors may be used to endorse or promote products derived
 *    from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

 /*
  * See: https://www.piday.org/million/
  */

#include "MonteCarloPi.h"
#include <algorithm>
#define CUDA_DRIVER_API
#include <helper_cuda.h>
#include <iostream>

#define ROUND_UP_TO_GRANULARITY(x, n) (((x + n - 1) / n) * n)

  // `ipcHandleTypeFlag` specifies the platform specific handle type this sample
  // uses for importing and exporting memory allocation. On Linux this sample
  // specifies the type as CU_MEM_HANDLE_TYPE_POSIX_FILE_DESCRIPTOR meaning that
  // file descriptors will be used. On Windows this sample specifies the type as
  // CU_MEM_HANDLE_TYPE_WIN32 meaning that NT HANDLEs will be used. The
  // ipcHandleTypeFlag variable is a convenience variable and is passed by value
  // to individual requests.
#if defined(__linux__)
CUmemAllocationHandleType ipcHandleTypeFlag = CU_MEM_HANDLE_TYPE_POSIX_FILE_DESCRIPTOR;
#else
CUmemAllocationHandleType ipcHandleTypeFlag = CU_MEM_HANDLE_TYPE_WIN32;
#endif

// Windows-specific LPSECURITYATTRIBUTES
void getDefaultSecurityDescriptor(CUmemAllocationProp *prop) {
#if defined(__linux__)
    return;
#elif defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64)
    static const char sddl[] = "D:P(OA;;GARCSDWDWOCCDCLCSWLODTWPRPCRFA;;;WD)";
    static OBJECT_ATTRIBUTES objAttributes;
    static bool objAttributesConfigured = false;

    if (!objAttributesConfigured) {
        PSECURITY_DESCRIPTOR secDesc;
        BOOL result = ConvertStringSecurityDescriptorToSecurityDescriptorA(
            sddl, SDDL_REVISION_1, &secDesc, NULL);
        if (result == 0) {
            printf("IPC failure: getDefaultSecurityDescriptor Failed! (%d)\n",
                GetLastError());
        }

        InitializeObjectAttributes(&objAttributes, NULL, 0, NULL, secDesc);

        objAttributesConfigured = true;
    }

    prop->win32HandleMetaData = &objAttributes;
    return;
#endif
}

__global__ void monte_carlo_kernel(vec2 *xyVector, float *pointsInsideCircle, float *numPointsInCircle, unsigned int numPoints, float time)
{
    const size_t stride = gridDim.x * blockDim.x;
    size_t tid = blockIdx.x * blockDim.x + threadIdx.x;
    float count = 0.0f;

    curandState rgnState;
    curand_init((unsigned long long)time, tid, 0, &rgnState);

    for (; tid < numPoints; tid += stride) {
        float x = curand_uniform(&rgnState);
        float y = curand_uniform(&rgnState);
        x = (2.0f * x) - 1.0f;
        y = (2.0f * y) - 1.0f;
        xyVector[tid][0] = x;
        xyVector[tid][1] = y;

        // Compute the distance of this point form the center(0, 0)
        float dist = sqrtf((x*x) + (y*y));

        // If distance is less than the radius of the unit circle, the point lies in the circle.
        pointsInsideCircle[tid] = (dist <= 1.0f);
        count += (dist <= 1.0f);
    }
    atomicAdd(numPointsInCircle, count);
}

MonteCarloPiSimulation::MonteCarloPiSimulation(size_t num_points) :
    m_xyVector(nullptr),
    m_pointsInsideCircle(nullptr),
    m_totalPointsInsideCircle(0),
    m_totalPointsSimulated(0),
    m_numPoints(num_points)
{
}

MonteCarloPiSimulation::~MonteCarloPiSimulation()
{
    if (m_numPointsInCircle) {
        checkCudaErrors(cudaFree(m_numPointsInCircle));
        m_numPointsInCircle = nullptr;
    }
    if (m_hostNumPointsInCircle) {
        checkCudaErrors(cudaFreeHost(m_hostNumPointsInCircle));
        m_hostNumPointsInCircle = nullptr;
    }

    cleanupSimulationAllocations();
}

void MonteCarloPiSimulation::initSimulation(int cudaDevice, cudaStream_t stream)
{
    m_cudaDevice = cudaDevice;
    getIdealExecutionConfiguration();

    // Allocate a position buffer that contains random location of the points in XY cartesian plane.
    // Allocate a bitmap buffer which holds information of whether a point in the position buffer is inside the unit circle or not.
    setupSimulationAllocations();

    checkCudaErrors(cudaMalloc((float **)&m_numPointsInCircle, sizeof(*m_numPointsInCircle)));
    checkCudaErrors(cudaMallocHost((float **)&m_hostNumPointsInCircle, sizeof(*m_hostNumPointsInCircle)));
}

void MonteCarloPiSimulation::stepSimulation(float time, cudaStream_t stream)
{

    checkCudaErrors(cudaMemsetAsync(m_numPointsInCircle, 0, sizeof(*m_numPointsInCircle), stream));

    monte_carlo_kernel << < m_blocks, m_threads, 0, stream >> > (m_xyVector, m_pointsInsideCircle, m_numPointsInCircle, m_numPoints, time);
    getLastCudaError("Failed to launch CUDA simulation");

    checkCudaErrors(cudaMemcpyAsync(m_hostNumPointsInCircle, m_numPointsInCircle, sizeof(*m_numPointsInCircle), cudaMemcpyDeviceToHost, stream));

    // Queue up a stream callback to compute and print the PI value.
    checkCudaErrors(cudaLaunchHostFunc(stream, this->computePiCallback, (void *)this));
}

void MonteCarloPiSimulation::computePiCallback(void *args)
{
    MonteCarloPiSimulation *cbData = (MonteCarloPiSimulation *)args;
    cbData->m_totalPointsInsideCircle += *(cbData->m_hostNumPointsInCircle);
    cbData->m_totalPointsSimulated += cbData->m_numPoints;
    double piValue = 4.0 * ((double)cbData->m_totalPointsInsideCircle / (double)cbData->m_totalPointsSimulated);
    printf("Approximate Pi value for %zd data points: %lf \n", cbData->m_totalPointsSimulated, piValue);
}

void MonteCarloPiSimulation::getIdealExecutionConfiguration()
{
    int warpSize = 0;
    int multiProcessorCount = 0;

    checkCudaErrors(cudaSetDevice(m_cudaDevice));
    checkCudaErrors(cudaDeviceGetAttribute(&warpSize, cudaDevAttrWarpSize, m_cudaDevice));

    // We don't need large block sizes, since there's not much inter-thread communication
    m_threads = warpSize;

    // Use the occupancy calculator and fill the gpu as best as we can
    checkCudaErrors(cudaOccupancyMaxActiveBlocksPerMultiprocessor(&m_blocks, monte_carlo_kernel, warpSize, 0));

    checkCudaErrors(cudaDeviceGetAttribute(&multiProcessorCount, cudaDevAttrMultiProcessorCount, m_cudaDevice));
    m_blocks *= multiProcessorCount;

    // Go ahead and the clamp the blocks to the minimum needed for this height/width
    m_blocks = std::min(m_blocks, (int)((m_numPoints + m_threads - 1) / m_threads));
}

void MonteCarloPiSimulation::setupSimulationAllocations()
{
    CUdeviceptr d_ptr = 0U;
    size_t granularity = 0;
    CUmemGenericAllocationHandle cudaPositionHandle, cudaInCircleHandle;

    CUmemAllocationProp allocProp = { };
    allocProp.type = CU_MEM_ALLOCATION_TYPE_PINNED;
    allocProp.location.type = CU_MEM_LOCATION_TYPE_DEVICE;
    allocProp.location.id = m_cudaDevice;
    allocProp.win32HandleMetaData = NULL;
    allocProp.requestedHandleTypes = ipcHandleTypeFlag;

    // Windows-specific LPSECURITYATTRIBUTES is required when
    // CU_MEM_HANDLE_TYPE_WIN32 is used. The security attribute defines the scope
    // of which exported allocations may be tranferred to other processes. For all
    // other handle types, pass NULL.
    getDefaultSecurityDescriptor(&allocProp);

    // Get the recommended granularity for m_cudaDevice.
    checkCudaErrors(cuMemGetAllocationGranularity(&granularity, &allocProp, CU_MEM_ALLOC_GRANULARITY_RECOMMENDED));

    size_t xyPositionVecSize = m_numPoints * sizeof(*m_xyVector);
    size_t inCircleVecSize = m_numPoints * sizeof(*m_pointsInsideCircle);

    size_t xyPositionSize = ROUND_UP_TO_GRANULARITY(xyPositionVecSize, granularity);
    size_t inCircleSize = ROUND_UP_TO_GRANULARITY(inCircleVecSize, granularity);
    m_totalAllocationSize = (xyPositionSize + inCircleSize);

    // Reserve the required contiguous VA space for the allocations
    checkCudaErrors(cuMemAddressReserve(&d_ptr, m_totalAllocationSize, granularity, 0U, 0));

    // Create the allocations as a pinned allocation on this device.
    // Create an allocation to store all the positions of points on the xy plane and a second
    // allocation which stores information if the corresponding position is inside the unit circle or not.
    checkCudaErrors(cuMemCreate(&cudaPositionHandle, xyPositionSize, &allocProp, 0));
    checkCudaErrors(cuMemCreate(&cudaInCircleHandle, inCircleSize, &allocProp, 0));

    // Export the allocation to a platform-specific handle. The type of handle
    // requested here must match the requestedHandleTypes field in the prop
    // structure passed to cuMemCreate. The handle obtained here will be passed to vulkan
    // to import the allocation.
    checkCudaErrors(cuMemExportToShareableHandle((void *)&m_posShareableHandle, cudaPositionHandle, ipcHandleTypeFlag, 0));
    checkCudaErrors(cuMemExportToShareableHandle((void *)&m_inCircleShareableHandle, cudaInCircleHandle, ipcHandleTypeFlag, 0));

    CUdeviceptr va_position = d_ptr;
    CUdeviceptr va_InCircle = va_position + xyPositionSize;
    m_pointsInsideCircle = (float *)va_InCircle;
    m_xyVector = (vec2 *)va_position;

    // Assign the chunk to the appropriate VA range
    checkCudaErrors(cuMemMap(va_position, xyPositionSize, 0, cudaPositionHandle, 0));
    checkCudaErrors(cuMemMap(va_InCircle, inCircleSize, 0, cudaInCircleHandle, 0));

    // Release the handles for the allocation. Since the allocation is currently mapped to a VA range
    // with a previous call to cuMemMap the actual freeing of memory allocation will happen on an eventual call to
    // cuMemUnmap. Thus the allocation will be kept live until it is unmapped.
    checkCudaErrors(cuMemRelease(cudaPositionHandle));
    checkCudaErrors(cuMemRelease(cudaInCircleHandle));

    CUmemAccessDesc accessDescriptor = {};
    accessDescriptor.location.id = m_cudaDevice;
    accessDescriptor.location.type = CU_MEM_LOCATION_TYPE_DEVICE;
    accessDescriptor.flags = CU_MEM_ACCESS_FLAGS_PROT_READWRITE;

    // Apply the access descriptor to the whole VA range. Essentially enables Read-Write access to the range.
    checkCudaErrors(cuMemSetAccess(d_ptr, m_totalAllocationSize, &accessDescriptor, 1));
}

void MonteCarloPiSimulation::cleanupSimulationAllocations()
{
    if (m_xyVector && m_pointsInsideCircle) {
        // Unmap the mapped virtual memory region
        // Since the handles to the mapped backing stores have already been released
        // by cuMemRelease, and these are the only/last mappings referencing them,
        // The backing stores will be freed.
        checkCudaErrors(cuMemUnmap((CUdeviceptr)m_xyVector, m_totalAllocationSize));

        checkIpcErrors(ipcCloseShareableHandle(m_posShareableHandle));
        checkIpcErrors(ipcCloseShareableHandle(m_inCircleShareableHandle));

        // Free the virtual address region.
        checkCudaErrors(cuMemAddressFree((CUdeviceptr)m_xyVector, m_totalAllocationSize));

        m_xyVector = nullptr;
        m_pointsInsideCircle = nullptr;
    }
}