#include

#include

#include

#include

#include

#include

#include

#define MATRIX_DIM 3

#define THRESHOLD 0.1

void LoadImage(const char *filename, float **data, int *height, int *width);

int SaveImage(const char *filename, float *image, int height, int width);

__host__ __device__

float GetValidPixelValue(float *image, int height, int width, int r, int c) {

// Prevent trying to read pixel locations outside the image

if (r < 0 || r >= height || c < 0 || c >= width)

return 0;

else

return image[r * width + c];

}

void SobelCPU(float *image, int height, int width, float *x_matrix,

float *y_matrix, float *output) {

for (int i = 0; i < height; i++) {

for (int j = 0; j < width; j++) {

float x_grad = 0;

x_grad += x_matrix[0 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i - 1, j - 1);

x_grad += x_matrix[0 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i - 1, j + 1);

x_grad += x_matrix[1 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i, j - 1);

x_grad += x_matrix[1 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i, j + 1);

x_grad += x_matrix[2 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i + 1, j - 1);

x_grad += x_matrix[2 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i + 1, j + 1);

float y_grad = 0;

y_grad += y_matrix[0 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i - 1, j - 1);

y_grad += y_matrix[0 * MATRIX_DIM + 1] *

GetValidPixelValue(image, height, width, i - 1, j);

y_grad += y_matrix[0 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i - 1, j + 1);

y_grad += y_matrix[2 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i + 1, j - 1);

y_grad += y_matrix[2 * MATRIX_DIM + 1] *

GetValidPixelValue(image, height, width, i + 1, j);

y_grad += y_matrix[2 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i + 1, j + 1);

float magnitude =

sqrt(x_grad / 8 * x_grad / 8 +

y_grad / 8 * y_grad / 8); // normalize gradients by dividing by 8

if (magnitude > 1) { // clamp to 1

output[i * width + j] = 1;

} else if (magnitude > THRESHOLD) {

output[i * width + j] = magnitude;

} else {

output[i * width + j] = 0;

}

}

}

}

__global__

void SobelGPUShared(float *image, int height, int width, float *x_matrix,

float *y_matrix, float *output) {

extern __shared__ float arr[];

int i = blockIdx.x*blockDim.x+threadIdx.x;

int j = blockIdx.y*blockDim.y+threadIdx.y;

//arr[i*j] = image[i*j];

if(i < height && j < width)

{

float x_grad=0;

x_grad += x_matrix[0 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i - 1, j - 1);

x_grad += x_matrix[0 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i - 1, j + 1);

x_grad += x_matrix[1 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i, j - 1);

x_grad += x_matrix[1 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i, j + 1);

x_grad += x_matrix[2 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i + 1, j - 1);

x_grad += x_matrix[2 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i + 1, j + 1);

__syncthreads();

float y_grad = 0;

y_grad += y_matrix[0 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i - 1, j - 1);

y_grad += y_matrix[0 * MATRIX_DIM + 1] *

GetValidPixelValue(image, height, width, i - 1, j);

y_grad += y_matrix[0 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i - 1, j + 1);

y_grad += y_matrix[2 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i + 1, j - 1);

y_grad += y_matrix[2 * MATRIX_DIM + 1] *

GetValidPixelValue(image, height, width, i + 1, j);

y_grad += y_matrix[2 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i + 1, j + 1);

__syncthreads();

float magnitude =

sqrt(x_grad / 8 * x_grad / 8 +

y_grad / 8 * y_grad / 8); // normalize gradients by dividing by 8

if (magnitude > 1) { // clamp to 1

output[i * width + j] = 1;

} else if (magnitude > THRESHOLD) {

output[i * width + j] = magnitude;

} else {

output[i * width + j] = 0;

}

// Implement this function, define a GPU kernel above it

}

}

__global__

void SobelGPU(float *image, int height, int width, float *x_matrix,

float *y_matrix, float *output) {

int i = blockIdx.x*blockDim.x+threadIdx.x;

int j = blockIdx.y*blockDim.y+threadIdx.y;

if(i < height && j < width)

{

float x_grad=0;

x_grad += x_matrix[0 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i - 1, j - 1);

x_grad += x_matrix[0 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i - 1, j + 1);

x_grad += x_matrix[1 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i, j - 1);

x_grad += x_matrix[1 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i, j + 1);

x_grad += x_matrix[2 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i + 1, j - 1);

x_grad += x_matrix[2 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i + 1, j + 1);

float y_grad = 0;

y_grad += y_matrix[0 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i - 1, j - 1);

y_grad += y_matrix[0 * MATRIX_DIM + 1] *

GetValidPixelValue(image, height, width, i - 1, j);

y_grad += y_matrix[0 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i - 1, j + 1);

y_grad += y_matrix[2 * MATRIX_DIM + 0] *

GetValidPixelValue(image, height, width, i + 1, j - 1);

y_grad += y_matrix[2 * MATRIX_DIM + 1] *

GetValidPixelValue(image, height, width, i + 1, j);

y_grad += y_matrix[2 * MATRIX_DIM + 2] *

GetValidPixelValue(image, height, width, i + 1, j + 1);

float magnitude =

sqrt(x_grad / 8 * x_grad / 8 +

y_grad / 8 * y_grad / 8); // normalize gradients by dividing by 8

if (magnitude > 1) { // clamp to 1

output[i * width + j] = 1;

} else if (magnitude > THRESHOLD) {

output[i * width + j] = magnitude;

} else {

output[i * width + j] = 0;

}

// Implement this function, define a GPU kernel above it

}

}

double WallTime() { // returns time in MS as a double

struct timeval tv;

gettimeofday(&tv, 0);

return tv.tv_sec * 1000.0 + tv.tv_usec / 1000.0;

}

void CheckResults(int height, int width, float *gpu_result, float *cpu_result){

for (int i = 0; i < height * width; i++) {

if (fabs(gpu_result[i] - cpu_result[i]) > 0.01) {

printf(" Mismatch at index %d: GPU: %f CPU %f ", i, gpu_result[i], cpu_result[i]);

return;

}

}

printf(" Ok!!! ");

}

int main(int argc, char *argv[]) {

float *data;

double start, end;

int height, width;

LoadImage("input.bmp", &data, &height, &width);

// CPU/GPU Data

float *x_matrix;

float *y_matrix;

float *image;

float *cpu_result;

float *gpu_result;

cudaMallocManaged(&image, sizeof(float) * height * width);

cudaMallocManaged(&x_matrix, sizeof(float) * MATRIX_DIM * MATRIX_DIM);

cudaMallocManaged(&y_matrix, sizeof(float) * MATRIX_DIM * MATRIX_DIM);

cudaMallocManaged(&cpu_result, sizeof(float) * height * width);

cudaMallocManaged(&gpu_result, sizeof(float) * height * width);

memcpy(image, data, sizeof(float) * height * width);

// Sobel x kernel matrix (already flipped)

x_matrix[0] = -1;

x_matrix[1] = 0;

x_matrix[2] = 1;

x_matrix[3] = -2;

x_matrix[4] = 0;

x_matrix[5] = 2;

x_matrix[6] = -1;

x_matrix[7] = 0;

x_matrix[7] = 1;

// Sobel y kernel matrix (already flipped)

y_matrix[0] = -1;

y_matrix[1] = -2;

y_matrix[2] = -1;

y_matrix[3] = 0;

y_matrix[4] = 0;

y_matrix[5] = 0;

y_matrix[6] = 1;

y_matrix[7] = 2;

y_matrix[8] = 1;

start = WallTime();

SobelCPU(image, height, width, x_matrix, y_matrix, cpu_result);

end = WallTime();

printf("CPU Runtime = %f msecs ", end - start);

dim3 threadsPerBlock(16,16);

dim3 numBlocks(height/threadsPerBlock.x, width/threadsPerBlock.y);

cudaMemcpy(image, data, sizeof(float)*height*width,cudaMemcpyHostToDevice);

SobelGPUShared<<>>(image, height, width, x_matrix, y_matrix, gpu_result);

cudaMemcpy(gpu_result,gpu_result, sizeof(float)*height*width,cudaMemcpyDeviceToHost);

cudaDeviceSynchronize();

CheckResults(height, width, gpu_result, cpu_result);

SaveImage("cpu_result.bmp", cpu_result, height, width);

SaveImage("gpu_result.bmp", gpu_result, height, width);

free(data);

cudaFree(image);

cudaFree(cpu_result);

cudaFree(gpu_result);

cudaFree(x_matrix);

cudaFree(y_matrix);

return 0;

}

Given this piece of Cuda code, how would you fix the shared memory sobel filter function SobelGPUShared so that the image is divided up into tiles that fit into shared memory?