I try to run the following code which would compute the dot product of two vectors, and the code can run well when the input number of GPU is 1, that is, the Omp package isn't really used, but when the number of GPU is 2, the GPU result is always 0, I don't know where is wrong, I just use usual parallel reduction in gpu code, and the seperate the work in N GPUs. I've check the code of multiGPUs run well when I don't use parallel reduction in gpu code, that is, I let C[i] = A[i]+B[i] and compute the sum at host.
// using multiple GPUs with OpenMP
// Includes
#include <stdio.h>
#include <stdlib.h>
#include <omp.h> // header for OpenMP
#include <cuda_runtime.h>
// Variables
float* h_A; // host vectors
float* h_B;
float* h_C;
float* h_D;
// Functions
void RandomInit(float*, int);
// Device code
__global__ void VecAdd(const float* A, const float* B, float* C, int N)
{
extern __shared__ float cache[];
int i = blockDim.x * blockIdx.x + threadIdx.x;
int cacheIndex = threadIdx.x;
float temp = 0.0; // register for each thread
while (i < N) {
temp += A[i]*B[i];
i += blockDim.x*gridDim.x;
}
cache[cacheIndex] = temp; // set the cache value
__syncthreads();
// perform parallel reduction, threadsPerBlock must be 2^m
int ib = blockDim.x/2;
while (ib != 0) {
if(cacheIndex < ib)
cache[cacheIndex] += cache[cacheIndex + ib];
__syncthreads();
ib /=2;
}
if(cacheIndex == 0)
C[blockIdx.x] = cache[0];
}
// Host code
int main(void)
{
printf("\n");
printf("Vector Dot Product with multiple GPUs \n");
int N, NGPU, cpu_thread_id=0;
int *Dev;
long mem = 1024*1024*1024; // 4 Giga for float data type.
printf("Enter the number of GPUs: ");
scanf("%d", &NGPU);
printf("%d\n", NGPU);
Dev = (int *)malloc(sizeof(int)*NGPU);
int numDev = 0;
printf("GPU device number: ");
for(int i = 0; i < NGPU; i++) {
scanf("%d", &Dev[i]);
printf("%d ",Dev[i]);
numDev++;
if(getchar() == '\n') break;
}
printf("\n");
if(numDev != NGPU) {
fprintf(stderr,"Should input %d GPU device numbers\n", NGPU);
exit(1);
}
printf("Enter the size of the vectors: ");
scanf("%d", &N);
printf("%d\n", N);
if (3*N > mem) {
printf("The size of these 3 vectors cannot be fitted into 4 Gbyte\n");
exit(1);
}
long size = N*sizeof(float);
// Set the sizes of threads and blocks
int threadsPerBlock;
printf("Enter the number of threads per block: ");
scanf("%d", &threadsPerBlock);
printf("%d\n", threadsPerBlock);
if(threadsPerBlock > 1024) {
printf("The number of threads per block must be less than 1024 ! \n");
exit(1);
}
int blocksPerGrid = (N + threadsPerBlock*NGPU - 1) / (threadsPerBlock*NGPU);
printf("The number of blocks is %d\n", blocksPerGrid);
if(blocksPerGrid > 2147483647) {
printf("The number of blocks must be less than 2147483647 ! \n");
exit(1);
}
long sb = blocksPerGrid*sizeof(float);
long sm = threadsPerBlock*sizeof(float);
// Allocate input vectors h_A and h_B in host memory
h_A = (float*)malloc(size);
h_B = (float*)malloc(size);
h_C = (float*)malloc(sb);
if (! h_A || ! h_B || ! h_C) {
printf("!!! Not enough memory.\n");
exit(1);
}
// Initialize input vectors
RandomInit(h_A, N);
RandomInit(h_B, N);
// declare cuda event for timer
cudaEvent_t start, stop;
// cudaEventCreate(&start); // events must be created after devices are set
// cudaEventCreate(&stop);
float Intime,gputime,Outime;
double h_G = 0.0;
omp_set_num_threads(NGPU);
#pragma omp parallel private(cpu_thread_id)
{
float *d_A, *d_B, *d_C;
cpu_thread_id = omp_get_thread_num();
cudaSetDevice(Dev[cpu_thread_id]);
// cudaSetDevice(cpu_thread_id);
// start the timer
if(cpu_thread_id == 0) {
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start,0);
}
// Allocate vectors in device memory
cudaMalloc((void**)&d_A, size/NGPU);
cudaMalloc((void**)&d_B, size/NGPU);
cudaMalloc((void**)&d_C, sb/NGPU);
// Copy vectors from host memory to device memory
cudaMemcpy(d_A, h_A+N/NGPU*cpu_thread_id, size/NGPU, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B+N/NGPU*cpu_thread_id, size/NGPU, cudaMemcpyHostToDevice);
#pragma omp barrier
// stop the timer
if(cpu_thread_id == 0) {
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
cudaEventElapsedTime( &Intime, start, stop);
printf("Data input time for GPU: %f (ms) \n",Intime);
}
// start the timer
if(cpu_thread_id == 0) cudaEventRecord(start,0);
VecAdd<<<blocksPerGrid, threadsPerBlock, sm>>>(d_A, d_B, d_C, N/NGPU);
cudaDeviceSynchronize();
// stop the timer
if(cpu_thread_id == 0) {
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
cudaEventElapsedTime( &gputime, start, stop);
printf("Processing time for GPU: %f (ms) \n",gputime);
printf("GPU Gflops: %f\n",3*N/(1000000.0*gputime));
}
// Copy result from device memory to host memory
// h_C contains the result in host memory
// start the timer
if(cpu_thread_id == 0) cudaEventRecord(start,0);
cudaMemcpy(h_C+blocksPerGrid/NGPU*cpu_thread_id, d_C, sb/NGPU, cudaMemcpyDeviceToHost);
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
//compute the solution
for (int i = 0; i < blocksPerGrid; i++) {
h_G += (double) h_C[i];
}
// stop the timer
if(cpu_thread_id == 0) {
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
cudaEventElapsedTime( &Outime, start, stop);
printf("Data output time for GPU: %f (ms) \n",Outime);
}
}
float gputime_tot;
gputime_tot = Intime + gputime + Outime;
printf("Total time for GPU: %f (ms) \n",gputime_tot);
// start the timer
cudaEventRecord(start,0);
double h_D = 0.0; // compute the reference solution
for (int i = 0; i < N; ++i)
h_D += (double) h_A[i]*h_B[i];
// h_D[i] = 1.0/cos(h_A[i]) + 1.0/cos(h_B[i]);
// stop the timer
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
float cputime;
cudaEventElapsedTime( &cputime, start, stop);
printf("Processing time for CPU: %f (ms) \n",cputime);
printf("CPU Gflops: %f\n",3*N/(1000000.0*cputime));
printf("Speed up of GPU = %f\n", cputime/gputime_tot);
// Destroy timer
cudaEventDestroy(start);
cudaEventDestroy(stop);
// check result
printf("Check result:\n");
// for (int i = 0; i < N; ++i) {
// diff = abs(h_D[i] - h_C[i]);
// sum += diff*diff;
// }
double diff = abs( (h_D - h_G)/h_D );
printf("|(h_G - h_D)/h_D|=%20.15e\n",diff);
printf("h_G =%20.15e\n",h_G);
printf("h_D =%20.15e\n",h_D);
for (int i=0; i < NGPU; i++) {
cudaSetDevice(i);
cudaDeviceReset();
}
return 0;
}
// Allocates an array with random float entries.
void RandomInit(float* data, int n)
{
for (int i = 0; i < n; ++i)
data[i] = rand() / (float)RAND_MAX;
}
First, its good practice to use proper CUDA error checking.
Clearly the work needs to be divided by the number of GPUs. But its unclear what your variables should mean. Let's drive a stake in the ground and say that blocksPerGrid will be the definition of the number of blocks in the kernel launch (for each GPU). That's consistent (at least) with your actual kernel invocations as you have shown them.
If we start there, then blocksPerGrid is going to be "multiplied" (i.e. scaled up) by the number of GPUs in order to cover your entire problem size. Let's go through your code and "harmonize" the calculations. For example, for two GPUs, a vector size of 1048576, and 512 threads per block, we expect blocksPerGrid to be 1024, because 2x1024x512 = 1048576. This is consistent with your calculation of blocksPerGrid itself and your kernel invocation.
1. This is incorrect:
long sb = blocksPerGrid*sizeof(float);
...
h_C = (float*)malloc(sb);
the host storage for the result needs to match (at least) the problem size. It needs to be one float item per block, times the number of GPUs. But sb is the storage size per GPU. We need to multiply it by the number of GPUs, when calculating the needed size for h_C.
2. This is incorrect:
cudaMalloc((void**)&d_C, sb/NGPU);
sb is already the storage size per GPU due to your calculation of blocksPerGrid. You should not divide it again by NGPU. When you do so, you now have threadblocks in each GPU that are attempting to write results to non-existent allocation, and your kernel would perform illegal behavior. Given a big enough problem and/or use of compute-sanitizer you would certainly witness this with the proper CUDA error checking I mentioned.
3. This is incorrect:
cudaMemcpy(h_C+blocksPerGrid/NGPU*cpu_thread_id, d_C, sb/NGPU, cudaMemcpyDeviceToHost);
for reasons we have already covered. h_C needs to cover the entire problem size, and the problem size per GPU is already covered by blocksPerGrid. It should not be further divided by NGPU, and sb is already the scaled per GPU, it should not be further divided by NGPU.
4. This is incorrect:
for (int i = 0; i < blocksPerGrid; i++) {
h_G += (double) h_C[i];
}
We have already covered the fact that your blocksPerGrid calculation is inherently a per-GPU calculation. It does not cover the whole problem size in the multi-GPU case.
5.
Your placement of the calculation of h_G is incorrect. We require that all OMP threads complete their work, before h_G result is calculated. Therefore this calculation needs to be after the closure of the OMP parallel region, to guarantee that all threads have updated their portion of h_C.
The following code has changes to address those issues, and seems to run correctly for me. To avoid user input and uncertainty, I have hard-coded some input values and changed the random initialization to one that is easy to assess for correctness:
$ cat t3.cu
#include <stdio.h>
#include <stdlib.h>
#include <omp.h> // header for OpenMP
#include <cuda_runtime.h>
// Variables
float* h_A; // host vectors
float* h_B;
float* h_C;
float* h_D;
// Functions
void RandomInit(float*, int);
// Device code
__global__ void VecAdd(const float* A, const float* B, float* C, int N)
{
extern __shared__ float cache[];
int i = blockDim.x * blockIdx.x + threadIdx.x;
int cacheIndex = threadIdx.x;
float temp = 0.0; // register for each thread
while (i < N) {
temp += A[i]*B[i];
i += blockDim.x*gridDim.x;
}
cache[cacheIndex] = temp; // set the cache value
__syncthreads();
// perform parallel reduction, threadsPerBlock must be 2^m
int ib = blockDim.x/2;
while (ib != 0) {
if(cacheIndex < ib)
cache[cacheIndex] += cache[cacheIndex + ib];
__syncthreads();
ib /=2;
}
if(cacheIndex == 0)
C[blockIdx.x] = cache[0];
}
// Host code
int main(void)
{
printf("\n");
printf("Vector Dot Product with multiple GPUs \n");
int N, NGPU, cpu_thread_id=0;
int *Dev;
long mem = 1024*1024*1024; // 4 Giga for float data type.
printf("Enter the number of GPUs: ");
//scanf("%d", &NGPU);
NGPU = 2;
printf("%d\n", NGPU);
Dev = (int *)malloc(sizeof(int)*NGPU);
int numDev = 0;
printf("GPU device number: ");
for(int i = 0; i < NGPU; i++) {
//scanf("%d", &Dev[i]);
Dev[i] = i;
printf("%d ",Dev[i]);
numDev++;
// if(getchar() == '\n') break;
}
printf("\n");
if(numDev != NGPU) {
fprintf(stderr,"Should input %d GPU device numbers\n", NGPU);
exit(1);
}
printf("Enter the size of the vectors: ");
//scanf("%d", &N);
N = 1048576;
printf("%d\n", N);
if (3*N > mem) {
printf("The size of these 3 vectors cannot be fitted into 4 Gbyte\n");
exit(1);
}
long size = N*sizeof(float);
// Set the sizes of threads and blocks
int threadsPerBlock;
printf("Enter the number of threads per block: ");
//scanf("%d", &threadsPerBlock);
threadsPerBlock = 512;
printf("%d\n", threadsPerBlock);
if(threadsPerBlock > 1024) {
printf("The number of threads per block must be less than 1024 ! \n");
exit(1);
}
int blocksPerGrid = (N + threadsPerBlock*NGPU - 1) / (threadsPerBlock*NGPU);
printf("The number of blocks is %d\n", blocksPerGrid);
if(blocksPerGrid > 2147483647) {
printf("The number of blocks must be less than 2147483647 ! \n");
exit(1);
}
long sb = blocksPerGrid*sizeof(float);
long sm = threadsPerBlock*sizeof(float);
// Allocate input vectors h_A and h_B in host memory
h_A = (float*)malloc(size);
h_B = (float*)malloc(size);
h_C = (float*)malloc(sb*NGPU);
if (! h_A || ! h_B || ! h_C) {
printf("!!! Not enough memory.\n");
exit(1);
}
// Initialize input vectors
RandomInit(h_A, N);
RandomInit(h_B, N);
// declare cuda event for timer
cudaEvent_t start, stop;
// cudaEventCreate(&start); // events must be created after devices are set
// cudaEventCreate(&stop);
float Intime,gputime,Outime;
double h_G = 0.0;
omp_set_num_threads(NGPU);
#pragma omp parallel private(cpu_thread_id)
{
float *d_A, *d_B, *d_C;
cpu_thread_id = omp_get_thread_num();
cudaSetDevice(Dev[cpu_thread_id]);
// cudaSetDevice(cpu_thread_id);
// start the timer
if(cpu_thread_id == 0) {
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start,0);
}
// Allocate vectors in device memory
cudaMalloc((void**)&d_A, size/NGPU);
cudaMalloc((void**)&d_B, size/NGPU);
cudaMalloc((void**)&d_C, sb);
// Copy vectors from host memory to device memory
cudaMemcpy(d_A, h_A+N/NGPU*cpu_thread_id, size/NGPU, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B+N/NGPU*cpu_thread_id, size/NGPU, cudaMemcpyHostToDevice);
#pragma omp barrier
// stop the timer
if(cpu_thread_id == 0) {
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
cudaEventElapsedTime( &Intime, start, stop);
printf("Data input time for GPU: %f (ms) \n",Intime);
}
// start the timer
if(cpu_thread_id == 0) cudaEventRecord(start,0);
VecAdd<<<blocksPerGrid, threadsPerBlock, sm>>>(d_A, d_B, d_C, N/NGPU);
cudaDeviceSynchronize();
// stop the timer
if(cpu_thread_id == 0) {
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
cudaEventElapsedTime( &gputime, start, stop);
printf("Processing time for GPU: %f (ms) \n",gputime);
printf("GPU Gflops: %f\n",3*N/(1000000.0*gputime));
}
// Copy result from device memory to host memory
// h_C contains the result in host memory
// start the timer
if(cpu_thread_id == 0) cudaEventRecord(start,0);
cudaMemcpy(h_C+blocksPerGrid*cpu_thread_id, d_C, sb, cudaMemcpyDeviceToHost);
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
// stop the timer
if(cpu_thread_id == 0) {
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
cudaEventElapsedTime( &Outime, start, stop);
printf("Data output time for GPU: %f (ms) \n",Outime);
}
}
//compute the solution
for (int i = 0; i < blocksPerGrid*NGPU; i++) {
h_G += (double) h_C[i];
}
float gputime_tot;
gputime_tot = Intime + gputime + Outime;
printf("Total time for GPU: %f (ms) \n",gputime_tot);
// start the timer
cudaEventRecord(start,0);
double h_D = 0.0; // compute the reference solution
for (int i = 0; i < N; ++i)
h_D += (double) h_A[i]*h_B[i];
// h_D[i] = 1.0/cos(h_A[i]) + 1.0/cos(h_B[i]);
// stop the timer
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
float cputime;
cudaEventElapsedTime( &cputime, start, stop);
printf("Processing time for CPU: %f (ms) \n",cputime);
printf("CPU Gflops: %f\n",3*N/(1000000.0*cputime));
printf("Speed up of GPU = %f\n", cputime/gputime_tot);
// Destroy timer
cudaEventDestroy(start);
cudaEventDestroy(stop);
// check result
printf("Check result:\n");
// for (int i = 0; i < N; ++i) {
// diff = abs(h_D[i] - h_C[i]);
// sum += diff*diff;
// }
double diff = abs( (h_D - h_G)/h_D );
printf("|(h_G - h_D)/h_D|=%20.15e\n",diff);
printf("h_G =%20.15e\n",h_G);
printf("h_D =%20.15e\n",h_D);
for (int i=0; i < NGPU; i++) {
cudaSetDevice(i);
cudaDeviceReset();
}
return 0;
}
// Allocates an array with random float entries.
void RandomInit(float* data, int n)
{
for (int i = 0; i < n; ++i)
data[i] = 1.0f; //rand() / (float)RAND_MAX;
}
$ nvcc -o t3 t3.cu -Xcompiler -fopenmp
$ compute-sanitizer ./t3
========= COMPUTE-SANITIZER
Vector Dot Product with multiple GPUs
Enter the number of GPUs: 2
GPU device number: 0 1
Enter the size of the vectors: 1048576
Enter the number of threads per block: 512
The number of blocks is 1024
Data input time for GPU: 2.405280 (ms)
Processing time for GPU: 8.202272 (ms)
GPU Gflops: 0.383519
Data output time for GPU: 0.429728 (ms)
Total time for GPU: 11.037280 (ms)
Processing time for CPU: 2.361696 (ms)
CPU Gflops: 1.331978
Speed up of GPU = 0.213974
Check result:
|(h_G - h_D)/h_D|=0.000000000000000e+00
h_G =1.048576000000000e+06
h_D =1.048576000000000e+06
========= ERROR SUMMARY: 0 errors
$
I'm not suggesting I have discovered every possible error in your code. The only test case I have tried is the one depicted.