reduction example for other than powers of 2

I have 2 arrays of 19600 elements. I am using a single kernel to multiply these arrays, do a sum reduction and also find the maximum.

I instantiate the kernel with 16 threads and 1225 blocks. Using the reduction example, I am able to do a sum reduction resulting in 1225 values. I am trying to understand how I can also sum 25 blocks such that I have 49 reduced values?

__shared__ float prod[threads];

int arrayIndex = gridDim.x*blockDim.x*blockIdx.y +blockDim.x*blockIdx.x + threadIdx.x;

c[arrayIndex] = a[arrayIndex]*b[arrayIndex];	

prod[threadIdx.x] = c[arrayIndex];

__syncthreads();

//reduce the values in each block //

if (threadIdx.x == 0)

{

	sumr[blockIdx.x] = prod[threadIdx.x];				

}

sumr is an array of 1225 values. I tried looping through the block Indexes, but that gives me the wrong value. Any suggestions?

You got at least two different working pieces of code for this in your other thread on the same subject. Why don’t you just use one of those?

You got at least two different working pieces of code for this in your other thread on the same subject. Why don’t you just use one of those?

@avidday: In my previous thread, I solved my problem, where I had to do reduction on all the values in an array efficiently. Here I am trying to do a reduce-by-key type where I want to only reduce each 25 values of the 1225 values obtained resulting in total of 49 values. So the key for reduction is the blockIdx.x and that’s Y I wrote this as a separate question.

I also thought abt ur suggestion in the previous thread regarding launching a kernel with the no. of blocks u want. However, if I am going to launch a second kernel to do this task, I would launch it with 49 blocks and 25 threads, But then since 25 is not a power of 2, I am not sure about how the reduction will work. Same as the scenario, when I launched the kernel with 49 blocks and 400 threads, the intermediate results stored in each block was wrong. Perhaps the reduction method relies so much on powers of 2, that I am lost when dealing with threads other than powers of 2.

@avidday: In my previous thread, I solved my problem, where I had to do reduction on all the values in an array efficiently. Here I am trying to do a reduce-by-key type where I want to only reduce each 25 values of the 1225 values obtained resulting in total of 49 values. So the key for reduction is the blockIdx.x and that’s Y I wrote this as a separate question.

I also thought abt ur suggestion in the previous thread regarding launching a kernel with the no. of blocks u want. However, if I am going to launch a second kernel to do this task, I would launch it with 49 blocks and 25 threads, But then since 25 is not a power of 2, I am not sure about how the reduction will work. Same as the scenario, when I launched the kernel with 49 blocks and 400 threads, the intermediate results stored in each block was wrong. Perhaps the reduction method relies so much on powers of 2, that I am lost when dealing with threads other than powers of 2.

In that case just use 32 threads per block, and zero the unused shared memory locations. Something like this:

__shared__ float vv80[32];

vv80[threadIdx.x] = (threadIdx.x < 25) ? input[blockstart+threadIdx.x] : 0.f;

// Normal shared memory reduction follows.....

In that case just use 32 threads per block, and zero the unused shared memory locations. Something like this:

__shared__ float vv80[32];

vv80[threadIdx.x] = (threadIdx.x < 25) ? input[blockstart+threadIdx.x] : 0.f;

// Normal shared memory reduction follows.....

I have an array of 20K values and I am reducing it over 50 blocks with 400 threads each.

My code looks like this:

getmax <<< num_blocks,block_size >>> (d_in, d_out1, d_indices);

__global__ void getmax(float *in1, float *out1, int *index)

{

	// Declare arrays to be in shared memory.

	__shared__ float max[threads];

	

	int nTotalThreads = blockDim.x;	// Total number of active threads

	float temp;

	float max_val;

	int max_index;

	int arrayIndex;

	

	// Calculate which element this thread reads from memory

	arrayIndex = gridDim.x*blockDim.x*blockIdx.y + blockDim.x*blockIdx.x + threadIdx.x;

	max[threadIdx.x] = in1[arrayIndex];

	max_val = max[threadIdx.x];

	max_index = blockDim.x*blockIdx.x + threadIdx.x;

	__syncthreads();

	while(nTotalThreads > 1)

	{

		int halfPoint = (nTotalThreads >> 1);

		if (threadIdx.x < halfPoint) 

		{

			temp = max[threadIdx.x + halfPoint];

			if (temp > max[threadIdx.x]) 

			{

				max[threadIdx.x] = temp;

				max_val = max[threadIdx.x];			

			}

		}

		__syncthreads();

		nTotalThreads = (nTotalThreads >> 1);	// divide by two.

	}

	if (threadIdx.x == 0)

	{

		out1[num_blocks*blockIdx.y + blockIdx.x] = max[threadIdx.x];

	}

	if(max[blockIdx.x] == max_val )

	{

		index[blockIdx.x] = max_index;	

	}

}

Th problem/issue here is that at some point “nTotalThreads” is not exactly a power of 2, resulting in garbage value for the index. The array out1 gives me the maximum value in each block, which is correct and validated. But the value of the index is wrong. For example: the max value in the first block occurs at index=40, but the kernel gives the values of index as 15. Similarly the value of the max in the second block is at 440, but the kernel gives 416.

Any suggestions??

@avidday: Don’t know hot to apply your suggestion regarding

vv80[threadIdx.x] = (threadIdx.x < 25) ? input[blockstart+threadIdx.x] : 0.f;

in this case?

I have an array of 20K values and I am reducing it over 50 blocks with 400 threads each.

My code looks like this:

getmax <<< num_blocks,block_size >>> (d_in, d_out1, d_indices);

__global__ void getmax(float *in1, float *out1, int *index)

{

	// Declare arrays to be in shared memory.

	__shared__ float max[threads];

	

	int nTotalThreads = blockDim.x;	// Total number of active threads

	float temp;

	float max_val;

	int max_index;

	int arrayIndex;

	

	// Calculate which element this thread reads from memory

	arrayIndex = gridDim.x*blockDim.x*blockIdx.y + blockDim.x*blockIdx.x + threadIdx.x;

	max[threadIdx.x] = in1[arrayIndex];

	max_val = max[threadIdx.x];

	max_index = blockDim.x*blockIdx.x + threadIdx.x;

	__syncthreads();

	while(nTotalThreads > 1)

	{

		int halfPoint = (nTotalThreads >> 1);

		if (threadIdx.x < halfPoint) 

		{

			temp = max[threadIdx.x + halfPoint];

			if (temp > max[threadIdx.x]) 

			{

				max[threadIdx.x] = temp;

				max_val = max[threadIdx.x];			

			}

		}

		__syncthreads();

		nTotalThreads = (nTotalThreads >> 1);	// divide by two.

	}

	if (threadIdx.x == 0)

	{

		out1[num_blocks*blockIdx.y + blockIdx.x] = max[threadIdx.x];

	}

	if(max[blockIdx.x] == max_val )

	{

		index[blockIdx.x] = max_index;	

	}

}

Th problem/issue here is that at some point “nTotalThreads” is not exactly a power of 2, resulting in garbage value for the index. The array out1 gives me the maximum value in each block, which is correct and validated. But the value of the index is wrong. For example: the max value in the first block occurs at index=40, but the kernel gives the values of index as 15. Similarly the value of the max in the second block is at 440, but the kernel gives 416.

Any suggestions??

@avidday: Don’t know hot to apply your suggestion regarding

vv80[threadIdx.x] = (threadIdx.x < 25) ? input[blockstart+threadIdx.x] : 0.f;

in this case?

suggestions/pointers?

suggestions/pointers?