/* Copyright (c) 2022, 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.
*/
/*
* This sample implements a conjugate gradient solver on GPU using
* Multi Block Cooperative Groups, also uses Unified Memory.
*
*/
// includes, system
#include <cuda_runtime.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// Utilities and system includes
#include <cooperative_groups.h>
#include <cooperative_groups/reduce.h>
#include <helper_cuda.h> // helper function CUDA error checking and initialization
#include <helper_functions.h> // helper for shared functions common to CUDA Samples
namespace cg = cooperative_groups;
const char *sSDKname = "conjugateGradientMultiBlockCG";
#define ENABLE_CPU_DEBUG_CODE 0
#define THREADS_PER_BLOCK 512
/* genTridiag: generate a random tridiagonal symmetric matrix */
void genTridiag(int *I, int *J, float *val, int N, int nz)
{
I[0] = 0, J[0] = 0, J[1] = 1;
val[0] = static_cast<float>(rand()) / RAND_MAX + 10.0f;
val[1] = static_cast<float>(rand()) / RAND_MAX;
int start;
for (int i = 1; i < N; i++) {
if (i > 1) {
I[i] = I[i - 1] + 3;
}
else {
I[1] = 2;
}
start = (i - 1) * 3 + 2;
J[start] = i - 1;
J[start + 1] = i;
if (i < N - 1) {
J[start + 2] = i + 1;
}
val[start] = val[start - 1];
val[start + 1] = static_cast<float>(rand()) / RAND_MAX + 10.0f;
if (i < N - 1) {
val[start + 2] = static_cast<float>(rand()) / RAND_MAX;
}
}
I[N] = nz;
}
// I - contains location of the given non-zero element in the row of the matrix
// J - contains location of the given non-zero element in the column of the
// matrix val - contains values of the given non-zero elements of the matrix
// inputVecX - input vector to be multiplied
// outputVecY - resultant vector
void cpuSpMV(int *I, int *J, float *val, int nnz, int num_rows, float alpha, float *inputVecX, float *outputVecY)
{
for (int i = 0; i < num_rows; i++) {
int num_elems_this_row = I[i + 1] - I[i];
float output = 0.0;
for (int j = 0; j < num_elems_this_row; j++) {
output += alpha * val[I[i] + j] * inputVecX[J[I[i] + j]];
}
outputVecY[i] = output;
}
return;
}
double dotProduct(float *vecA, float *vecB, int size)
{
double result = 0.0;
for (int i = 0; i < size; i++) {
result = result + (vecA[i] * vecB[i]);
}
return result;
}
void scaleVector(float *vec, float alpha, int size)
{
for (int i = 0; i < size; i++) {
vec[i] = alpha * vec[i];
}
}
void saxpy(float *x, float *y, float a, int size)
{
for (int i = 0; i < size; i++) {
y[i] = a * x[i] + y[i];
}
}
void cpuConjugateGrad(int *I, int *J, float *val, float *x, float *Ax, float *p, float *r, int nnz, int N, float tol)
{
int max_iter = 10000;
float alpha = 1.0;
float alpham1 = -1.0;
float r0 = 0.0, b, a, na;
cpuSpMV(I, J, val, nnz, N, alpha, x, Ax);
saxpy(Ax, r, alpham1, N);
float r1 = dotProduct(r, r, N);
int k = 1;
while (r1 > tol * tol && k <= max_iter) {
if (k > 1) {
b = r1 / r0;
scaleVector(p, b, N);
saxpy(r, p, alpha, N);
}
else {
for (int i = 0; i < N; i++)
p[i] = r[i];
}
cpuSpMV(I, J, val, nnz, N, alpha, p, Ax);
float dot = dotProduct(p, Ax, N);
a = r1 / dot;
saxpy(p, x, a, N);
na = -a;
saxpy(Ax, r, na, N);
r0 = r1;
r1 = dotProduct(r, r, N);
printf("\nCPU code iteration = %3d, residual = %e\n", k, sqrt(r1));
k++;
}
}
__device__ void gpuSpMV(int *I,
int *J,
float *val,
int nnz,
int num_rows,
float alpha,
float *inputVecX,
float *outputVecY,
cg::thread_block &cta,
const cg::grid_group &grid)
{
for (int i = grid.thread_rank(); i < num_rows; i += grid.size()) {
int row_elem = I[i];
int next_row_elem = I[i + 1];
int num_elems_this_row = next_row_elem - row_elem;
float output = 0.0;
for (int j = 0; j < num_elems_this_row; j++) {
// I or J or val arrays - can be put in shared memory
// as the access is random and reused in next calls of gpuSpMV function.
output += alpha * val[row_elem + j] * inputVecX[J[row_elem + j]];
}
outputVecY[i] = output;
}
}
__device__ void gpuSaxpy(float *x, float *y, float a, int size, const cg::grid_group &grid)
{
for (int i = grid.thread_rank(); i < size; i += grid.size()) {
y[i] = a * x[i] + y[i];
}
}
__device__ void gpuDotProduct(float *vecA,
float *vecB,
double *result,
int size,
const cg::thread_block &cta,
const cg::grid_group &grid)
{
extern __shared__ double tmp[];
double temp_sum = 0.0;
for (int i = grid.thread_rank(); i < size; i += grid.size()) {
temp_sum += static_cast<double>(vecA[i] * vecB[i]);
}
cg::thread_block_tile<32> tile32 = cg::tiled_partition<32>(cta);
temp_sum = cg::reduce(tile32, temp_sum, cg::plus<double>());
if (tile32.thread_rank() == 0) {
tmp[tile32.meta_group_rank()] = temp_sum;
}
cg::sync(cta);
if (tile32.meta_group_rank() == 0) {
temp_sum = tile32.thread_rank() < tile32.meta_group_size() ? tmp[tile32.thread_rank()] : 0.0;
temp_sum = cg::reduce(tile32, temp_sum, cg::plus<double>());
if (tile32.thread_rank() == 0) {
atomicAdd(result, temp_sum);
}
}
}
__device__ void gpuCopyVector(float *srcA, float *destB, int size, const cg::grid_group &grid)
{
for (int i = grid.thread_rank(); i < size; i += grid.size()) {
destB[i] = srcA[i];
}
}
__device__ void
gpuScaleVectorAndSaxpy(const float *x, float *y, float a, float scale, int size, const cg::grid_group &grid)
{
for (int i = grid.thread_rank(); i < size; i += grid.size()) {
y[i] = a * x[i] + scale * y[i];
}
}
extern "C" __global__ void gpuConjugateGradient(int *I,
int *J,
float *val,
float *x,
float *Ax,
float *p,
float *r,
double *dot_result,
int nnz,
int N,
float tol)
{
cg::thread_block cta = cg::this_thread_block();
cg::grid_group grid = cg::this_grid();
int max_iter = 10000;
float alpha = 1.0;
float alpham1 = -1.0;
float r0 = 0.0, r1, b, a, na;
gpuSpMV(I, J, val, nnz, N, alpha, x, Ax, cta, grid);
cg::sync(grid);
gpuSaxpy(Ax, r, alpham1, N, grid);
cg::sync(grid);
gpuDotProduct(r, r, dot_result, N, cta, grid);
cg::sync(grid);
r1 = *dot_result;
int k = 1;
while (r1 > tol * tol && k <= max_iter) {
if (k > 1) {
b = r1 / r0;
gpuScaleVectorAndSaxpy(r, p, alpha, b, N, grid);
}
else {
gpuCopyVector(r, p, N, grid);
}
cg::sync(grid);
gpuSpMV(I, J, val, nnz, N, alpha, p, Ax, cta, grid);
if (threadIdx.x == 0 && blockIdx.x == 0)
*dot_result = 0.0;
cg::sync(grid);
gpuDotProduct(p, Ax, dot_result, N, cta, grid);
cg::sync(grid);
a = r1 / *dot_result;
gpuSaxpy(p, x, a, N, grid);
na = -a;
gpuSaxpy(Ax, r, na, N, grid);
r0 = r1;
cg::sync(grid);
if (threadIdx.x == 0 && blockIdx.x == 0)
*dot_result = 0.0;
cg::sync(grid);
gpuDotProduct(r, r, dot_result, N, cta, grid);
cg::sync(grid);
r1 = *dot_result;
k++;
}
}
bool areAlmostEqual(float a, float b, float maxRelDiff)
{
float diff = fabsf(a - b);
float abs_a = fabsf(a);
float abs_b = fabsf(b);
float largest = abs_a > abs_b ? abs_a : abs_b;
if (diff <= largest * maxRelDiff) {
return true;
}
else {
printf("maxRelDiff = %.8e\n", maxRelDiff);
printf("diff %.8e > largest * maxRelDiff %.8e therefore %.8e and %.8e are not "
"same\n",
diff,
largest * maxRelDiff,
a,
b);
return false;
}
}
int main(int argc, char **argv)
{
int N = 0, nz = 0, *I = NULL, *J = NULL;
float *val = NULL;
const float tol = 1e-5f;
float *x;
float *rhs;
float r1;
float *r, *p, *Ax;
cudaEvent_t start, stop;
printf("Starting [%s]...\n", sSDKname);
// This will pick the best possible CUDA capable device
cudaDeviceProp deviceProp;
int devID = findCudaDevice(argc, (const char **)argv);
checkCudaErrors(cudaGetDeviceProperties(&deviceProp, devID));
if (!deviceProp.managedMemory) {
// This sample requires being run on a device that supports Unified Memory
fprintf(stderr, "Unified Memory not supported on this device\n");
exit(EXIT_WAIVED);
}
// This sample requires being run on a device that supports Cooperative Kernel
// Launch
if (!deviceProp.cooperativeLaunch) {
printf("\nSelected GPU (%d) does not support Cooperative Kernel Launch, "
"Waiving the run\n",
devID);
exit(EXIT_WAIVED);
}
// Statistics about the GPU device
printf("> GPU device has %d Multi-Processors, SM %d.%d compute capabilities\n\n",
deviceProp.multiProcessorCount,
deviceProp.major,
deviceProp.minor);
/* Generate a random tridiagonal symmetric matrix in CSR format */
N = 1048576;
nz = (N - 2) * 3 + 4;
cudaMallocManaged(reinterpret_cast<void **>(&I), sizeof(int) * (N + 1));
cudaMallocManaged(reinterpret_cast<void **>(&J), sizeof(int) * nz);
cudaMallocManaged(reinterpret_cast<void **>(&val), sizeof(float) * nz);
genTridiag(I, J, val, N, nz);
cudaMallocManaged(reinterpret_cast<void **>(&x), sizeof(float) * N);
cudaMallocManaged(reinterpret_cast<void **>(&rhs), sizeof(float) * N);
double *dot_result;
cudaMallocManaged(reinterpret_cast<void **>(&dot_result), sizeof(double));
*dot_result = 0.0;
// temp memory for CG
checkCudaErrors(cudaMallocManaged(reinterpret_cast<void **>(&r), N * sizeof(float)));
checkCudaErrors(cudaMallocManaged(reinterpret_cast<void **>(&p), N * sizeof(float)));
checkCudaErrors(cudaMallocManaged(reinterpret_cast<void **>(&Ax), N * sizeof(float)));
cudaDeviceSynchronize();
checkCudaErrors(cudaEventCreate(&start));
checkCudaErrors(cudaEventCreate(&stop));
#if ENABLE_CPU_DEBUG_CODE
float *Ax_cpu = reinterpret_cast<float *>(malloc(sizeof(float) * N));
float *r_cpu = reinterpret_cast<float *>(malloc(sizeof(float) * N));
float *p_cpu = reinterpret_cast<float *>(malloc(sizeof(float) * N));
float *x_cpu = reinterpret_cast<float *>(malloc(sizeof(float) * N));
for (int i = 0; i < N; i++) {
r_cpu[i] = 1.0;
Ax_cpu[i] = x_cpu[i] = 0.0;
}
#endif
for (int i = 0; i < N; i++) {
r[i] = rhs[i] = 1.0;
x[i] = 0.0;
}
void *kernelArgs[] = {
(void *)&I,
(void *)&J,
(void *)&val,
(void *)&x,
(void *)&Ax,
(void *)&p,
(void *)&r,
(void *)&dot_result,
(void *)&nz,
(void *)&N,
(void *)&tol,
};
int sMemSize = sizeof(double) * ((THREADS_PER_BLOCK / 32) + 1);
int numBlocksPerSm = 0;
int numThreads = THREADS_PER_BLOCK;
checkCudaErrors(
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&numBlocksPerSm, gpuConjugateGradient, numThreads, sMemSize));
int numSms = deviceProp.multiProcessorCount;
dim3 dimGrid(numSms * numBlocksPerSm, 1, 1), dimBlock(THREADS_PER_BLOCK, 1, 1);
checkCudaErrors(cudaEventRecord(start, 0));
checkCudaErrors(
cudaLaunchCooperativeKernel((void *)gpuConjugateGradient, dimGrid, dimBlock, kernelArgs, sMemSize, NULL));
checkCudaErrors(cudaEventRecord(stop, 0));
checkCudaErrors(cudaDeviceSynchronize());
float time;
checkCudaErrors(cudaEventElapsedTime(&time, start, stop));
r1 = *dot_result;
printf("GPU Final, residual = %e, kernel execution time = %f ms\n", sqrt(r1), time);
#if ENABLE_CPU_DEBUG_CODE
cpuConjugateGrad(I, J, val, x_cpu, Ax_cpu, p_cpu, r_cpu, nz, N, tol);
#endif
float rsum, diff, err = 0.0;
for (int i = 0; i < N; i++) {
rsum = 0.0;
for (int j = I[i]; j < I[i + 1]; j++) {
rsum += val[j] * x[J[j]];
}
diff = fabs(rsum - rhs[i]);
if (diff > err) {
err = diff;
}
}
checkCudaErrors(cudaFree(I));
checkCudaErrors(cudaFree(J));
checkCudaErrors(cudaFree(val));
checkCudaErrors(cudaFree(x));
checkCudaErrors(cudaFree(rhs));
checkCudaErrors(cudaFree(r));
checkCudaErrors(cudaFree(p));
checkCudaErrors(cudaFree(Ax));
checkCudaErrors(cudaFree(dot_result));
checkCudaErrors(cudaEventDestroy(start));
checkCudaErrors(cudaEventDestroy(stop));
#if ENABLE_CPU_DEBUG_CODE
free(Ax_cpu);
free(r_cpu);
free(p_cpu);
free(x_cpu);
#endif
printf("Test Summary: Error amount = %f \n", err);
fprintf(stdout, "&&&& conjugateGradientMultiBlockCG %s\n", (sqrt(r1) < tol) ? "PASSED" : "FAILED");
exit((sqrt(r1) < tol) ? EXIT_SUCCESS : EXIT_FAILURE);
}