I have a number (say a million) of small matrices 4 x 3. I would like to do several simple operations with them and I would like my CUDA kernel to parallelize the matrices index only, (not the row/column operations). Let me explain better: I pass as an input to my GPU Kernel an array of matrices A[MatrixNumb][row][col] and I would like operation parallelization to be only on the MatrixNumb (therefore I want to force the operation in one dimension. The example below is with 3 Matrices only, for simplicity. It compiles and runs, however it gives me the wrong results. Basically, it returns the same matrices I give it as an input. I cannot understand why and if I am making any mistake, how can I re-write/think the code? I wrote the code using also CudaMallocManaged, in order to have shared memory between host and device, however it gives me the same results using the classic CudaMalloc and using memcpy.
Source.cpp
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <iostream>
#include <assert.h>
#include <chrono>
#include <random>
#include <time.h>
#include <math.h>
#include <cuda_runtime.h>
#include "device_launch_parameters.h"
#include <cuda.h>
#include <device_functions.h>
using namespace std;
__global__ void SVD(double*** a, const int m, const int n, const int numMatrices, double** w)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
// I would like that each thread runs these loops independently
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
a[idx][i][j] = (a[idx][i][j] * a[idx][i][j]) * 3.14;
}
}
for (int j = 0; j < n; j++) {
w[idx][j] = 3.14 * a[idx][1][j]* a[idx][1][j];
}
}
int main()
{
const int n = 3;
const int m = 4;
const int lda = m;
const int numMatrices = 3;
random_device device;
mt19937 generator(device());
uniform_real_distribution<double> distribution(1., 5.);
// create pointers
double*** A = new double** [numMatrices];
double** w = new double* [numMatrices];
//ALLOCATE SHARED MEMORY
for (int nm = 0; nm < numMatrices; nm++) {
A[nm] = new double* [lda];
w[nm] = new double[n];
for (int i = 0; i < lda; i++) {
A[nm][i] = new double[n];
for (int j = 0; j < n; j++) {
cudaMallocManaged((void**)&A[nm][i][j], sizeof(double));
cudaMallocManaged((void**)&w[nm][j], sizeof(double));
}
}
}
cout << " memory allocated" << endl;
//FILL MATRICES INTO SHARED MEMORY
for (int nm = 0; nm < numMatrices; nm++) {
A[nm] = new double* [lda];
w[nm] = new double[n];
for (int i = 0; i < lda; i++) {
A[nm][i] = new double[n];
for (int j = 0; j < n; j++) {
A[nm][i][j] = distribution(generator);
w[nm][j] = 0.0;
}
}
}
cout << " matrix filled " << endl;
// PRINT MATRICES BEFORE CUDA OPERATION
for (int nm = 0; nm < numMatrices; nm++) {
for (int i = 0; i < lda; i++) {
for (int j = 0; j < n; j++) {
cout << A[nm][i][j] << " ";
}
cout << endl;
}
cout << endl;
}
//KERNEL ----------------------------------------------------------------------
int NThreads = 3;
int NBlocks = int(numMatrices / NThreads + 1);
SVD << <NBlocks, NThreads >> > (A, n, m, numMatrices, w);
cudaDeviceSynchronize();
cout << " Kernel done " << endl << endl;
cout << " --- GPU --- " << endl;
cout << " NEW MATRIX: " << endl;
for (int nm = 0; nm < numMatrices; nm++) {
for (int i = 0; i < lda; i++) {
for (int j = 0; j < n; j++) {
cout << A[nm][i][j] << " ";
}
cout << endl;
}
cout << endl;
}
cout << " NEW VECTOR RESULTS: " << endl;
for (int nm = 0; nm < numMatrices; nm++) {
for (int i = 0; i < n; i++) {
cout << w[nm][i] << " ";
}
cout << endl;
}
cout << endl;
//FREE THE DEVICE'S MEMORY -----------------------------------------------------
cudaFree(A);
cudaFree(w);
cout << " Cuda free " << endl << endl;
return 0;
}
The (wrong) output I get is the following:
memory allocated
matrix filled
1.28689 3.76588 3.88649
1.52547 4.42371 2.62566
1.48002 3.33719 1.58413
3.78243 2.8394 3.0249
1.14322 1.70261 2.02784
2.86852 2.87918 3.2896
4.87268 3.52447 1.58414
3.52306 3.84931 3.18212
1.76397 1.41317 4.9765
1.63338 4.79316 2.64009
1.99873 1.72617 1.15974
1.18922 4.21513 1.6695
Kernel done
--- GPU ---
NEW MATRIX:
1.28689 3.76588 3.88649
1.52547 4.42371 2.62566
1.48002 3.33719 1.58413
3.78243 2.8394 3.0249
1.14322 1.70261 2.02784
2.86852 2.87918 3.2896
4.87268 3.52447 1.58414
3.52306 3.84931 3.18212
1.76397 1.41317 4.9765
1.63338 4.79316 2.64009
1.99873 1.72617 1.15974
1.18922 4.21513 1.6695
NEW VECTOR RESULTS:
0 0 0
0 0 0
0 0 0
Cuda free
I expected to se the new matrices and the vectors modified by the operations of: a[idx][i][j] = (a[idx][i][j] * a[idx][i][j]) * 3.14; however, It looks like the code does not see the kernel or the kernel does not work properly.