aaron/ispc
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

added workable .cu files for stencil & mandelbrot

Browse Source
This commit is contained in:
Evghenii
2013-11-08 10:00:49 +01:00
parent cb7cbec0d5
commit 426afc7377
10 changed files with 645 additions and 160 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
examples_cuda/stencil/a.out
View File

Binary file not shown.

4
examples_cuda/stencil/stencil.cu
View File

@@ -1,6 +1,6 @@
#define programCount 32
#define programIndex threadIdx.x
#define taskIndex blockIdx.x
#define programIndex (threadIdx.x & 31)
#define taskIndex (blockIdx.x*4 + (threadIdx.x >> 5))
__device__ static void
stencil_step( int x0, int x1,

BIN
examples_cuda/stencil/stencil.cubin
View File

Binary file not shown.

66
examples_cuda/stencil/stencil.ispc
View File

@@ -34,13 +34,14 @@
#ifdef __NVPTX__
#warning "emitting DEVICE code"
#define taskIndex blockIndex0()
#define taskCount blockCount0()
#define programIndex laneIndex()
#define programCount warpSize()
#else
#warning "emitting HOST code"
#endif
static void
static inline void
stencil_step(uniform int x0, uniform int x1,
uniform int y0, uniform int y1,
uniform int z0, uniform int z1,
@@ -50,29 +51,62 @@ stencil_step(uniform int x0, uniform int x1,
const uniform int Nxy = Nx * Ny;
// foreach (z = z0 ... z1, y = y0 ... y1, x = x0 ... x1)
#if 0
#define VER1
#endif
#ifdef VER1
const uniform long x1o = 1;
const uniform long x2o = 2;
const uniform long x3o = 3;
const uniform long y1o = Nx;
const uniform long y2o = Nx*2;
const uniform long y3o = Nx*3;
const uniform long z1o = Nxy;
const uniform long z2o = Nxy*2;
const uniform long z3o = Nxy*3;
#endif
for (uniform int z = z0; z < z1; z++)
for (uniform int y = y0; y < y1; y++)
for (uniform int xb = x0; xb < x1; xb += programCount)
{
const int x = xb + programIndex;
int index = (z * Nxy) + (y * Nx) + x;
const int index_base = (z * Nxy) + (y * Nx);
for (uniform int xb = x0; xb < x1; xb += programCount)
{
const int x = xb + programIndex;
int index = index_base + x;
#ifndef VER1
#define A_cur(x, y, z) Ain[index + (x) + ((y) * Nx) + ((z) * Nxy)]
#define A_next(x, y, z) Aout[index + (x) + ((y) * Nx) + ((z) * Nxy)]
double div = coef[0] * A_cur(0, 0, 0) +
double div = coef[0] * A_cur(0, 0, 0) +
coef[1] * (A_cur(+1, 0, 0) + A_cur(-1, 0, 0) +
A_cur(0, +1, 0) + A_cur(0, -1, 0) +
A_cur(0, 0, +1) + A_cur(0, 0, -1)) +
A_cur(0, +1, 0) + A_cur(0, -1, 0) +
A_cur(0, 0, +1) + A_cur(0, 0, -1)) +
coef[2] * (A_cur(+2, 0, 0) + A_cur(-2, 0, 0) +
A_cur(0, +2, 0) + A_cur(0, -2, 0) +
A_cur(0, 0, +2) + A_cur(0, 0, -2)) +
A_cur(0, +2, 0) + A_cur(0, -2, 0) +
A_cur(0, 0, +2) + A_cur(0, 0, -2)) +
coef[3] * (A_cur(+3, 0, 0) + A_cur(-3, 0, 0) +
A_cur(0, +3, 0) + A_cur(0, -3, 0) +
A_cur(0, 0, +3) + A_cur(0, 0, -3));
A_cur(0, +3, 0) + A_cur(0, -3, 0) +
A_cur(0, 0, +3) + A_cur(0, 0, -3));
#else
#define A_cur(x, y, z) Ain [index + (x) + (y) + (z)]
#define A_next(x, y, z) Aout[index + (x) + (y) + (z)]
double div = coef[0] * A_cur(0, 0, 0) +
coef[1] * (A_cur(+x1o, 0, 0) + A_cur(-x1o, 0, 0) +
A_cur(0, +y1o, 0) + A_cur(0, -y1o, 0) +
A_cur(0, 0, +z1o) + A_cur(0, 0, -z1o)) +
coef[2] * (A_cur(+x2o, 0, 0) + A_cur(-x2o, 0, 0) +
A_cur(0, +y2o, 0) + A_cur(0, -y2o, 0) +
A_cur(0, 0, +z2o) + A_cur(0, 0, -z2o)) +
coef[3] * (A_cur(+x3o, 0, 0) + A_cur(-x3o, 0, 0) +
A_cur(0, +y3o, 0) + A_cur(0, -y3o, 0) +
A_cur(0, 0, +z3o) + A_cur(0, 0, -z3o));
#endif
if (x < x1)
A_next(0, 0, 0) = 2.0d0 * A_cur(0, 0, 0) - A_next(0, 0, 0) +
vsq[index] * div;
}
if (x < x1)
A_next(0, 0, 0) = 2.0d0 * A_cur(0, 0, 0) - A_next(0, 0, 0) +
vsq[index] * div;
}
}
}
@@ -83,6 +117,8 @@ stencil_step_task(uniform int x0, uniform int x1,
uniform int Nx, uniform int Ny, uniform int Nz,
uniform const double coef[4], uniform const double vsq[],
uniform const double Ain[], uniform double Aout[]) {
if(taskIndex >= taskCount) return;
stencil_step(x0, x1, y0, y1, z0+taskIndex, z0+taskIndex+1,
Nx, Ny, Nz, coef, vsq, Ain, Aout);
}

261
examples_cuda/stencil/stencil_cu.cpp
View File

@@ -132,11 +132,12 @@ void memcpyH2D(CUdeviceptr d_buf, void * h_buf, const size_t size)
checkCudaErrors(cuMemcpyHtoD(d_buf, h_buf, size));
}
#define deviceLaunch(func,nbx,nby,nbz,params) \
checkCudaErrors(cuFuncSetCacheConfig((func), CU_FUNC_CACHE_PREFER_L1)); \
checkCudaErrors( \
cuLaunchKernel( \
(func), \
(nbx), (nby), (nbz), \
32, 1, 1, \
((nbx-1)/(128/32)+1), (nby), (nbz), \
128, 1, 1, \
0, NULL, (params), NULL \
));
@@ -144,6 +145,38 @@ typedef CUdeviceptr devicePtr;
/**************/
#include <vector>
std::vector<char> readBinary(const char * filename)
{
std::vector<char> buffer;
FILE *fp = fopen(filename, "rb");
if (!fp )
{
fprintf(stderr, "file %s not found\n", filename);
assert(0);
}
#if 0
char c;
while ((c = fgetc(fp)) != EOF)
buffer.push_back(c);
#else
fseek(fp, 0, SEEK_END);
const unsigned long long size = ftell(fp); /*calc the size needed*/
fseek(fp, 0, SEEK_SET);
buffer.resize(size);
if (fp == NULL){ /*ERROR detection if file == empty*/
fprintf(stderr, "Error: There was an Error reading the file %s \n",filename);
exit(1);
}
else if (fread(&buffer[0], sizeof(char), size, fp) != size){ /* if count of read bytes != calculated size of .bin file -> ERROR*/
fprintf(stderr, "Error: There was an Error reading the file %s \n", filename);
exit(1);
}
#endif
fprintf(stderr, " read buffer of size= %d bytes \n", (int)buffer.size());
return buffer;
}
extern "C"
{
@@ -155,15 +188,21 @@ extern "C"
void CUDALaunch(
void **handlePtr,
const char * module_name,
const char * module,
const char * module_1,
const char * func_name,
void **func_args,
int countx, int county, int countz)
{
assert(module_name != NULL);
assert(module != NULL);
assert(module_1 != NULL);
assert(func_name != NULL);
assert(func_args != NULL);
#if 1
const char * module = module_1;
#else
const std::vector<char> module_str = readBinary("kernel.cubin");
const char * module = &module_str[0];
#endif
CUmodule cudaModule = loadModule(module);
CUfunction cudaFunction = getFunction(cudaModule, func_name);
deviceLaunch(cudaFunction, countx, county, countz, func_args);
@@ -184,134 +223,134 @@ extern "C"
extern void loop_stencil_serial(int t0, int t1, int x0, int x1,
int y0, int y1, int z0, int z1,
int Nx, int Ny, int Nz,
const double coef[5],
const double vsq[],
double Aeven[], double Aodd[]);
int y0, int y1, int z0, int z1,
int Nx, int Ny, int Nz,
const double coef[5],
const double vsq[],
double Aeven[], double Aodd[]);
void InitData(int Nx, int Ny, int Nz, double *A[2], double *vsq) {
int offset = 0;
for (int z = 0; z < Nz; ++z)
for (int y = 0; y < Ny; ++y)
for (int x = 0; x < Nx; ++x, ++offset) {
A[0][offset] = (x < Nx / 2) ? x / double(Nx) : y / double(Ny);
A[1][offset] = 0;
vsq[offset] = x*y*z / double(Nx * Ny * Nz);
}
int offset = 0;
for (int z = 0; z < Nz; ++z)
for (int y = 0; y < Ny; ++y)
for (int x = 0; x < Nx; ++x, ++offset) {
A[0][offset] = (x < Nx / 2) ? x / double(Nx) : y / double(Ny);
A[1][offset] = 0;
vsq[offset] = x*y*z / double(Nx * Ny * Nz);
}
}
int main() {
int Nx = 256, Ny = 256, Nz = 256;
int width = 4;
double *Aserial[2], *Aispc[2];
Aserial[0] = new double [Nx * Ny * Nz];
Aserial[1] = new double [Nx * Ny * Nz];
Aispc[0] = new double [Nx * Ny * Nz];
Aispc[1] = new double [Nx * Ny * Nz];
double *vsq = new double [Nx * Ny * Nz];
int Nx = 256, Ny = 256, Nz = 256;
int width = 4;
double *Aserial[2], *Aispc[2];
Aserial[0] = new double [Nx * Ny * Nz];
Aserial[1] = new double [Nx * Ny * Nz];
Aispc[0] = new double [Nx * Ny * Nz];
Aispc[1] = new double [Nx * Ny * Nz];
double *vsq = new double [Nx * Ny * Nz];
double coeff[4] = { 0.5, -.25, .125, -.0625 };
/*******************/
createContext();
/*******************/
double coeff[4] = { 0.5, -.25, .125, -.0625 };
const size_t bufsize = sizeof(double)*Nx*Ny*Nz;
devicePtr d_Aispc0 = deviceMalloc(bufsize);
devicePtr d_Aispc1 = deviceMalloc(bufsize);
devicePtr d_vsq = deviceMalloc(bufsize);
devicePtr d_coeff = deviceMalloc(4*sizeof(double));
/*******************/
createContext();
/*******************/
const size_t bufsize = sizeof(double)*Nx*Ny*Nz;
devicePtr d_Aispc0 = deviceMalloc(bufsize);
devicePtr d_Aispc1 = deviceMalloc(bufsize);
devicePtr d_vsq = deviceMalloc(bufsize);
devicePtr d_coeff = deviceMalloc(4*sizeof(double));
InitData(Nx, Ny, Nz, Aispc, vsq);
InitData(Nx, Ny, Nz, Aispc, vsq);
//
// Compute the image using the ispc implementation on one core; report
// the minimum time of three runs.
//
double minTimeISPC = 1e30;
for (int i = 0; i < 3; ++i) {
reset_and_start_timer();
loop_stencil_ispc(0, 6, width, Nx - width, width, Ny - width,
width, Nz - width, Nx, Ny, Nz, coeff, vsq,
Aispc[0], Aispc[1]);
double dt = get_elapsed_mcycles();
minTimeISPC = std::min(minTimeISPC, dt);
}
//
// Compute the image using the ispc implementation on one core; report
// the minimum time of three runs.
//
double minTimeISPC = 1e30;
for (int i = 0; i < 3; ++i) {
reset_and_start_timer();
loop_stencil_ispc(0, 6, width, Nx - width, width, Ny - width,
width, Nz - width, Nx, Ny, Nz, coeff, vsq,
Aispc[0], Aispc[1]);
double dt = get_elapsed_mcycles();
minTimeISPC = std::min(minTimeISPC, dt);
}
printf("[stencil ispc 1 core]:\t\t[%.3f] million cycles\n", minTimeISPC);
InitData(Nx, Ny, Nz, Aispc, vsq);
printf("[stencil ispc 1 core]:\t\t[%.3f] million cycles\n", minTimeISPC);
memcpyH2D(d_Aispc0, Aispc[0], bufsize);
memcpyH2D(d_Aispc1, Aispc[1], bufsize);
memcpyH2D(d_vsq, vsq, bufsize);
memcpyH2D(d_coeff, coeff, 4*sizeof(double));
//
// Compute the image using the ispc implementation with tasks; report
// the minimum time of three runs.
//
double minTimeISPCTasks = 1e30;
for (int i = 0; i < 3; ++i) {
reset_and_start_timer();
loop_stencil_ispc_tasks(0, 6, width, Nx - width, width, Ny - width,
width, Nz - width, Nx, Ny, Nz, (double*)d_coeff, (double*)d_vsq,
(double*)d_Aispc0, (double*)d_Aispc1);
double dt = get_elapsed_mcycles();
minTimeISPCTasks = std::min(minTimeISPCTasks, dt);
}
memcpyD2H(Aispc[1], d_Aispc1, bufsize);
//memcpyD2H(Aispc[1], d_vsq, bufsize);
InitData(Nx, Ny, Nz, Aispc, vsq);
printf("[stencil ispc + tasks]:\t\t[%.3f] million cycles\n", minTimeISPCTasks);
memcpyH2D(d_Aispc0, Aispc[0], bufsize);
memcpyH2D(d_Aispc1, Aispc[1], bufsize);
memcpyH2D(d_vsq, vsq, bufsize);
memcpyH2D(d_coeff, coeff, 4*sizeof(double));
//
// Compute the image using the ispc implementation with tasks; report
// the minimum time of three runs.
//
double minTimeISPCTasks = 1e30;
for (int i = 0; i < 3; ++i) {
reset_and_start_timer();
loop_stencil_ispc_tasks(0, 6, width, Nx - width, width, Ny - width,
width, Nz - width, Nx, Ny, Nz, (double*)d_coeff, (double*)d_vsq,
(double*)d_Aispc0, (double*)d_Aispc1);
double dt = get_elapsed_mcycles();
minTimeISPCTasks = std::min(minTimeISPCTasks, dt);
}
memcpyD2H(Aispc[1], d_Aispc1, bufsize);
//memcpyD2H(Aispc[1], d_vsq, bufsize);
InitData(Nx, Ny, Nz, Aserial, vsq);
printf("[stencil ispc + tasks]:\t\t[%.3f] million cycles\n", minTimeISPCTasks);
//
// And run the serial implementation 3 times, again reporting the
// minimum time.
//
double minTimeSerial = 1e30;
for (int i = 0; i < 3; ++i) {
reset_and_start_timer();
loop_stencil_serial(0, 6, width, Nx-width, width, Ny - width,
width, Nz - width, Nx, Ny, Nz, coeff, vsq,
Aserial[0], Aserial[1]);
double dt = get_elapsed_mcycles();
minTimeSerial = std::min(minTimeSerial, dt);
}
InitData(Nx, Ny, Nz, Aserial, vsq);
printf("[stencil serial]:\t\t[%.3f] million cycles\n", minTimeSerial);
//
// And run the serial implementation 3 times, again reporting the
// minimum time.
//
double minTimeSerial = 1e30;
for (int i = 0; i < 3; ++i) {
reset_and_start_timer();
loop_stencil_serial(0, 6, width, Nx-width, width, Ny - width,
width, Nz - width, Nx, Ny, Nz, coeff, vsq,
Aserial[0], Aserial[1]);
double dt = get_elapsed_mcycles();
minTimeSerial = std::min(minTimeSerial, dt);
}
printf("\t\t\t\t(%.2fx speedup from ISPC, %.2fx speedup from ISPC + tasks)\n",
minTimeSerial / minTimeISPC, minTimeSerial / minTimeISPCTasks);
printf("[stencil serial]:\t\t[%.3f] million cycles\n", minTimeSerial);
// Check for agreement
int offset = 0;
int nerr = 0;
for (int z = 0; z < Nz; ++z)
for (int y = 0; y < Ny; ++y)
for (int x = 0; x < Nx; ++x, ++offset) {
printf("\t\t\t\t(%.2fx speedup from ISPC, %.2fx speedup from ISPC + tasks)\n",
minTimeSerial / minTimeISPC, minTimeSerial / minTimeISPCTasks);
double error = fabsf((Aserial[1][offset] - Aispc[1][offset]) /
Aserial[1][offset]);
if (error > 1e-3)
{
if (nerr < 100)
printf("Error @ (%d,%d,%d): ispc = %g, serial = %g error= %g\n",
x, y, z, Aispc[1][offset], Aserial[1][offset], error);
nerr++;
}
}
// Check for agreement
int offset = 0;
int nerr = 0;
for (int z = 0; z < Nz; ++z)
for (int y = 0; y < Ny; ++y)
for (int x = 0; x < Nx; ++x, ++offset) {
fprintf(stderr, " nerr= %d frac= %g \n", nerr, 1.0*nerr/(1.0*Nx*Ny*Nz));
/*******************/
destroyContext();
/*******************/
double error = fabsf((Aserial[1][offset] - Aispc[1][offset]) /
Aserial[1][offset]);
if (error > 1e-3)
{
if (nerr < 100)
printf("Error @ (%d,%d,%d): ispc = %g, serial = %g error= %g\n",
x, y, z, Aispc[1][offset], Aserial[1][offset], error);
nerr++;
}
}
return 0;
fprintf(stderr, " nerr= %d frac= %g \n", nerr, 1.0*nerr/(1.0*Nx*Ny*Nz));
/*******************/
destroyContext();
/*******************/
return 0;
}