Paralel Reduction With less than 8000 values

[avidday]

For this code, the number of blocks must always be equal to the number of columns. This model is 1 block processes 1 column. If you have 6000 columns, you need to launch 6000 blocks. If you foresee using more than 65535 columns, then you will need to modify the kernel to use a 2D grid, which would allow you to have 4294836225 columns. If you ever had more columns than that, you will have to split the matrix up and process it is sections. That would also be the way to parallelise the problem over multiple GPUs.

[David_Lisin]

For this code, the number of blocks must always be equal to the number of columns. This model is 1 block processes 1 column. If you have 6000 columns, you need to launch 6000 blocks. If you foresee using more than 65535 columns, then you will need to modify the kernel to use a 2D grid, which would allow you to have 4294836225 columns. If you ever had more columns than that, you will have to split the matrix up and process it is sections. That would also be the way to parallelise the problem over multiple GPUs.

The matrix I work with (6000x65) consists fo 6000 rows and 65 columns, this i execute with:

kernelFunction<<<numColumns, 128, 128*sizeof(double)>>>(numRows,numRows,inMatrix,outMatrix);

But as said, results are slightly different…

I have no intention of using more than max numBlocks columns, so that is not a problem :)

If there is no other reason for slightly different results, I shall leave as is, although the change does seem strange, considering that secuential sums do give the correct result.

Thank you, David

[Sarnath]

Just going only by your last post, I think it is similar with a problem that I faced with Reduction:
[url=“The Official NVIDIA Forums | NVIDIA”]The Official NVIDIA Forums | NVIDIA

You may check it out and feel satisfied that one another person too has wondered what is wrong…

[avidday]

All I can suggest about the differences in results is to try coding a small test outside of matlab using normalized inputs. You should get within 1 last place digit, I would have thought Representation limitations should not come into play for 6000 numbers between 0 and 1, and the code is straight addition, where multiply-add fusion/intermediate result rounding cannot influence the results.

It is your project, your verification process, and your choice. Just for the record, I hereby declare that code was provided on an as-is basis, with no warranty, guarantees, or expressions of fitness for purpose, implied or otherwise. You use it at your own risk.

[ilovei7]

Hello avidday,

I am sorry. I still could not get your kernel running properly. I could not understand “lda” and “m”. What do they mean?

Would you please provide a short host code to run the kernel?

I am very very new on CUDA… so very very sorry for wasting your time here. But i really like to implement your efficient way to replace my naive code.

Many many Thanks!

I don’t think you understand what I am suggesting, or how the parallel reduction can be used in this case.

This is (I am guessing) your “naïve” approach, which you launch with 1 block:

__global__ void columnsumnaive(const int m, const int lda, 

				double *in, double *out)

{

	unsigned int idx = __mul24(threadIdx.x, lda);

	double sum	   = in[idx];

	// Loop through the column and sum

	// all the entries. (ssslllloooowwww.....)

	for(int i=1; i<m; i++)

		sum += in[idx+i];

	out[threadIdx.x] = sum;

}

and this is a shared memory parallel summation, with in-kernel reduction in shared memory which should be launched with one block per column:

__global__ void columnsumreduction(const int m, const int lda, 

				double *in, double *out)

{

	extern __shared__  double buff[];

	unsigned int tidx   = threadIdx.x;

	unsigned int idx	= tidx + __mul24(blockIdx.x, lda);

	unsigned int idxmax = m + __mul24(blockIdx.x, lda);

	buff[tidx] = 0.;

	// Compute partial sums down the column until we

	// have covered the whole column

	while (idx < idxmax) {

		buff[tidx] += in[idx];

		idx += blockDim.x;

	}

	__syncthreads();

	// Parallel reduction of the partial sum 

	// in shared memory

	if (blockDim.x == 512) {

		if (tidx < 256) 

			buff[tidx] += buff[tidx + 256];

		__syncthreads(); 

	}

	if (blockDim.x >= 256) {

		if (tidx < 128) 

			buff[tidx] += buff[tidx + 128];

		__syncthreads(); 

	}

	if (blockDim.x >= 128) {

		if (tidx < 64) 

			buff[tidx] += buff[tidx + 64];

		__syncthreads(); 

	}

	

	if (tidx < 32) {

		if (blockDim.x >=  64) {

			buff[tidx] += buff[tidx + 32];

		}

		if (blockDim.x >=  32) {

			buff[tidx] += buff[tidx + 16];

		}

		if (blockDim.x >=  16) {

			buff[tidx] += buff[tidx + 8];

		}

		if (blockDim.x >=   8) {

			buff[tidx] += buff[tidx + 4];

		}

		if (blockDim.x >=   4) {

			buff[tidx] += buff[tidx + 2];

		}

		if (blockDim.x >=   2) {

			buff[tidx] += buff[tidx + 1];

		}

	}

	

	// write result for this block to global mem 

	if (tidx == 0)

		out[blockIdx.x] = buff[0];

	

}

the difference in execution time is between the two of them is more than 40x on my GTX275 (one block with 65 threads for the naïve case, and 65 blocks with 128 threads for the reduction):

avid@cuda:~/code/columnsum$ ./columnsum 

Column summation on a random double precision array

For m=6000 n=65

Reduction summation: time=6.681600e-02 ms, deviation from cpu=0.000000e+00

Naive summation:	 time=2.820704e+00 ms, deviation from cpu=0.000000e+00

That is it takes less than 0.07ms to do your 6000 row x 65 column example in double precision for column-major ordered storage (and I should add there are probably pretty terrible shared memory bank conflicts and serialization in the reduction code I posted).

[avidday]

I could not understand “lda” and “m”. What do they mean?

lda is the leading dimension (or pitch) of the array and m is the column length.

[HickoryDock]

Hello avidday,

I am sorry. I still could not get your kernel running properly. I could not understand “lda” and “m”. What do they mean?

Would you please provide a short host code to run the kernel?

I am very very new on CUDA… so very very sorry for wasting your time here. But i really like to implement your efficient way to replace my naive code.

Many many Thanks!

Yes, it would be great to get some working code.

I compiled and ran, but just got a load of zeros… ;*(

I called it using:

void sumCols(Matrix C,Matrix A)

{

	int numBlocks=A.w;

		int numThreads=512;

		const int m=A.h*A.w;

		const int lda=A.h;

	   columnsumreduction<<<numThreads,numBlocks>>>(m,lda,A.data,C.data);

}

[HickoryDock]

Yes, it would be great to get some working code.

I compiled and ran, but just got a load of zeros… ;*(

I called it using:

void sumCols(Matrix C,Matrix A)

{

	int numBlocks=A.w;

		int numThreads=512;

		const int m=A.h*A.w;

		const int lda=A.h;

	   columnsumreduction<<<numThreads,numBlocks>>>(m,lda,A.data,C.data);

}

Schoolboy error: numThreads and numBlocks around the wrong way… Cheers!