I went ahead and wrote my own loop unrolling code, and attached it to
this message. Is it a crime against python?
(Note that this code isn't finished; I'm just looking for stylistic
feedback on the unrolling bit)
On Fri, Jun 12, 2009 at 5:26 PM, Andreas
Klöckner<[email protected]> wrote:
> On Freitag 12 Juni 2009, Andrew Wagner wrote:
>> Hello-
> SourceModule supports a no_extern_c keyword argument. However, once you use
> that, all the names in the CUDA module become "mangled" [1]. If you can live
> with just a few entrypoints that you manually declare extern "C", then this is
> likely a good way.
>
> [1] http://en.wikipedia.org/wiki/Name_mangling
Thanks! I'll try that if I run into some C++isms that are harder to
work around.
>> 2) If I just rip out the templates and do the loop unrolling in
>> python, are there some pythonic examples of doing the loop unrolling
>> out there? (I could roll my own, but I'm a novice so it might be ugly)
>
> http://documen.tician.de/pycuda/metaprog.html
>
> http://is.gd/10clQ
Thanks Andreas; I'll look more closely at that codepy stuff!
Cheers,
Drew
import numpy
#from helper_functions import *
#from plotting import *
import pycuda.autoinit
import pycuda.driver as cuda
import time
import string
KERNEL_RADIUS = 8
KERNEL_W = 2 * KERNEL_RADIUS + 1
template = '''
//24-bit multiplication is faster on G80,
//but we must be sure to multiply integers
//only within [-8M, 8M - 1] range
#define IMUL(a, b) __mul24(a, b)
////////////////////////////////////////////////////////////////////////////////
// Kernel configuration
////////////////////////////////////////////////////////////////////////////////
#define KERNEL_RADIUS $KERNEL_RADIUS
#define KERNEL_W $KERNEL_W
__device__ __constant__ float d_Kernel[KERNEL_W];
// Assuming ROW_TILE_W, KERNEL_RADIUS_ALIGNED and dataW
// are multiples of coalescing granularity size,
// all global memory operations are coalesced in convolutionRowGPU()
#define ROW_TILE_W 128
#define KERNEL_RADIUS_ALIGNED 16
// Assuming COLUMN_TILE_W and dataW are multiples
// of coalescing granularity size, all global memory operations
// are coalesced in convolutionColumnGPU()
#define COLUMN_TILE_W 16
#define COLUMN_TILE_H 48'''
# Ignore the ugly templated unrolling code...
'''
////////////////////////////////////////////////////////////////////////////////
// Loop unrolling templates, needed for best performance
////////////////////////////////////////////////////////////////////////////////
template<int i> __device__ float convolutionRow(float *data){
return
data[KERNEL_RADIUS - i] * d_Kernel[i]
+ convolutionRow<i - 1>(data);
}
template<> __device__ float convolutionRow<-1>(float *data){
return 0;
}
template<int i> __device__ float convolutionColumn(float *data){
return
data[(KERNEL_RADIUS - i) * COLUMN_TILE_W] * d_Kernel[i]
+ convolutionColumn<i - 1>(data);
}
template<> __device__ float convolutionColumn<-1>(float *data){
return 0;
}'''
template += '''
////////////////////////////////////////////////////////////////////////////////
// Row convolution filter
////////////////////////////////////////////////////////////////////////////////
__global__ void convolutionRowGPU(
float *d_Result,
float *d_Data,
int dataW,
int dataH
){
//Data cache
__shared__ float data[KERNEL_RADIUS + ROW_TILE_W + KERNEL_RADIUS];
//Current tile and apron limits, relative to row start
const int tileStart = IMUL(blockIdx.x, ROW_TILE_W);
const int tileEnd = tileStart + ROW_TILE_W - 1;
const int apronStart = tileStart - KERNEL_RADIUS;
const int apronEnd = tileEnd + KERNEL_RADIUS;
//Clamp tile and apron limits by image borders
const int tileEndClamped = min(tileEnd, dataW - 1);
const int apronStartClamped = max(apronStart, 0);
const int apronEndClamped = min(apronEnd, dataW - 1);
//Row start index in d_Data[]
const int rowStart = IMUL(blockIdx.y, dataW);
//Aligned apron start. Assuming dataW and ROW_TILE_W are multiples
//of half-warp size, rowStart + apronStartAligned is also a
//multiple of half-warp size, thus having proper alignment
//for coalesced d_Data[] read.
const int apronStartAligned = tileStart - KERNEL_RADIUS_ALIGNED;
const int loadPos = apronStartAligned + threadIdx.x;
//Set the entire data cache contents
//Load global memory values, if indices are within the image borders,
//or initialize with zeroes otherwise
if(loadPos >= apronStart){
const int smemPos = loadPos - apronStart;
data[smemPos] =
((loadPos >= apronStartClamped) && (loadPos <= apronEndClamped)) ?
d_Data[rowStart + loadPos] : 0;
}
//Ensure the completness of the loading stage
//because results, emitted by each thread depend on the data,
//loaded by another threads
__syncthreads();
const int writePos = tileStart + threadIdx.x;
//Assuming dataW and ROW_TILE_W are multiples of half-warp size,
//rowStart + tileStart is also a multiple of half-warp size,
//thus having proper alignment for coalesced d_Result[] write.
if(writePos <= tileEndClamped){
const int smemPos = writePos - apronStart;
float sum = 0;
'''
# Ignore ugly, broken loop unrolling
'''
#ifdef UNROLL_INNER
sum = convolutionRow<2 * KERNEL_RADIUS>(data + smemPos);
#else
'''
originalLoop = '''
for(int k = -KERNEL_RADIUS; k <= KERNEL_RADIUS; k++)
sum += data[smemPos + k] * d_Kernel[KERNEL_RADIUS - k];
'''
unrolledLoop = ''
for k in range(-KERNEL_RADIUS, KERNEL_RADIUS+1):
loopTemplate = string.Template('sum += data[smemPos + $k] * d_Kernel[KERNEL_RADIUS - $k];\n')
unrolledLoop += loopTemplate.substitute(k=k)
# A slightly less pythonic version?
#k = str(k)
#unrolledLoop += 'sum += data[smemPos + ' + k + '] * d_Kernel[KERNEL_RADIUS - ' + k + '];\n'
print unrolledLoop
template += originalLoop
template += unrolledLoop
template += '''
d_Result[rowStart + writePos] = sum;
}
}
////////////////////////////////////////////////////////////////////////////////
// Column convolution filter
////////////////////////////////////////////////////////////////////////////////
__global__ void convolutionColumnGPU(
float *d_Result,
float *d_Data,
int dataW,
int dataH,
int smemStride,
int gmemStride
){
//Data cache
__shared__ float data[COLUMN_TILE_W * (KERNEL_RADIUS + COLUMN_TILE_H + KERNEL_RADIUS)];
//Current tile and apron limits, in rows
const int tileStart = IMUL(blockIdx.y, COLUMN_TILE_H);
const int tileEnd = tileStart + COLUMN_TILE_H - 1;
const int apronStart = tileStart - KERNEL_RADIUS;
const int apronEnd = tileEnd + KERNEL_RADIUS;
//Clamp tile and apron limits by image borders
const int tileEndClamped = min(tileEnd, dataH - 1);
const int apronStartClamped = max(apronStart, 0);
const int apronEndClamped = min(apronEnd, dataH - 1);
//Current column index
const int columnStart = IMUL(blockIdx.x, COLUMN_TILE_W) + threadIdx.x;
//Shared and global memory indices for current column
int smemPos = IMUL(threadIdx.y, COLUMN_TILE_W) + threadIdx.x;
int gmemPos = IMUL(apronStart + threadIdx.y, dataW) + columnStart;
//Cycle through the entire data cache
//Load global memory values, if indices are within the image borders,
//or initialize with zero otherwise
for(int y = apronStart + threadIdx.y; y <= apronEnd; y += blockDim.y){
data[smemPos] =
((y >= apronStartClamped) && (y <= apronEndClamped)) ?
d_Data[gmemPos] : 0;
smemPos += smemStride;
gmemPos += gmemStride;
}
//Ensure the completness of the loading stage
//because results, emitted by each thread depend on the data,
//loaded by another threads
__syncthreads();
//Shared and global memory indices for current column
smemPos = IMUL(threadIdx.y + KERNEL_RADIUS, COLUMN_TILE_W) + threadIdx.x;
gmemPos = IMUL(tileStart + threadIdx.y , dataW) + columnStart;
//Cycle through the tile body, clamped by image borders
//Calculate and output the results
for(int y = tileStart + threadIdx.y; y <= tileEndClamped; y += blockDim.y){
float sum = 0;
#ifdef UNROLL_INNER
sum = convolutionColumn<2 * KERNEL_RADIUS>(data + smemPos);
#else
for(int k = -KERNEL_RADIUS; k <= KERNEL_RADIUS; k++)
sum +=
data[smemPos + IMUL(k, COLUMN_TILE_W)] *
d_Kernel[KERNEL_RADIUS - k];
#endif
d_Result[gmemPos] = sum;
smemPos += smemStride;
gmemPos += gmemStride;
}
}
'''
nvidia_separable_convolution_kernel_template = string.Template(template)
nvidia_separable_convolution_kernel_code = nvidia_separable_convolution_kernel_template.substitute(KERNEL_RADIUS = KERNEL_RADIUS, KERNEL_W = KERNEL_W)
nvidia_separable_convolution_kernel_module = cuda.SourceModule(nvidia_separable_convolution_kernel_code)
nvidia_separable_convolution_kernel_function_row = nvidia_separable_convolution_kernel_module.get_function('convolutionRowGPU')
nvidia_separable_convolution_kernel_function_col = nvidia_separable_convolution_kernel_module.get_function('convolutionColumnGPU')
# Helper functions for computing alignment...
def iDivUp(a, b):
# Round a / b to nearest higher integer value
a = numpy.int32(a)
b = numpy.int32(b)
return (a / b + 1) if (a % b != 0) else (a / b)
def iDivDown(a, b):
# Round a / b to nearest lower integer value
a = numpy.int32(a)
b = numpy.int32(b)
return a / b;
def iAlignUp(a, b):
# Align a to nearest higher multiple of b
a = numpy.int32(a)
b = numpy.int32(b)
return (a - a % b + b) if (a % b != 0) else a
def iAlignDown(a, b):
# Align a to nearest lower multiple of b
a = numpy.int32(a)
b = numpy.int32(b)
return a - a % b
def convolution_cuda(sourceImage, destImage, filterx, filtery):
# Perform separable convolution on sourceImage using CUDA.
DATA_W = iAlignUp(3072, 16);
DATA_H = 3072;
BYTES_PER_WORD = 4; # 4 for float32
DATA_SIZE = DATA_W * DATA_H * BYTES_PER_WORD;
KERNEL_SIZE = KERNEL_W * BYTES_PER_WORD;
# I'm not done here...
def test_convolution_cuda():
# Test the convolution kernel.
# I'm not done here either!...
pass
if __name__ == '__main__':
test_convolution_cuda()
boo = raw_input('Pausing so you can look at results... <Enter> to finish...')
_______________________________________________
PyCUDA mailing list
[email protected]
http://tiker.net/mailman/listinfo/pycuda_tiker.net
