#define PROGRAM_FILE_PATH "regression_analysis.cl"

#include <stdio.h>
#include <CL/cl.h>

/* CL Globals */
cl_device_id device;
cl_context context;
cl_program program;
cl_command_queue queue;
cl_mem int_buf_1, int_buf_2, int_buf_3, int_buf_4, d_mat_1, d_mat_2, d_mat_3, d_mat_4, d_mat_5, d_mat_6,
	d_mat_ptr_1, d_mat_ptr_2, d_vec_1, return_flt_buf_1, return_flt_buf_2, return_flt_buf_3,
	return_flt_buf_4, return_int_ptr_1, return_int_ptr_2, return_vec_1, return_vec_2,
	return_vec_3, return_flt_ptr_1, return_flt_ptr_2, d_vec_1, d_vec_2, d_vec_3, return_temp_buf;
cl_kernel kernel;

/* Find a GPU or CPU associated with the first available platform */
cl_device_id create_device() {

	cl_platform_id platform;
	cl_device_id dev;
	char name[128];
	int err;

	/* Identify a platform */
	err = clGetPlatformIDs(1, &platform, NULL);
	if (err < 0) {
		perror("Couldn't identify a platform");
		exit(1);
	}

	/* Access a device */
	err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_ACCELERATOR, 1, &dev, NULL);
	if (err == CL_DEVICE_NOT_FOUND) {
		err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &dev, NULL);
		if (err == CL_DEVICE_NOT_FOUND) {
			err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_CPU, 1, &dev, NULL);
			if (err < 0) {
				perror("Couldn't access any devices");
				exit(1);
			} else {
				clGetDeviceInfo(dev, CL_DEVICE_NAME, 128, name, NULL);
				fprintf(stdout, "Created a dispatch queue using the %s\n", name);
				printf("Accessing CPU... \n");
			}
		} else {
			clGetDeviceInfo(dev, CL_DEVICE_NAME, 128, name, NULL);
			fprintf(stdout, "Created a dispatch queue using the %s\n", name);
			printf("Accessing GPU... \n");
		}
	} else {
		clGetDeviceInfo(dev, CL_DEVICE_NAME, 128, name, NULL);
		fprintf(stdout, "Created a dispatch queue using the %s\n", name);
		printf("Accessing Accelerator... \n");
	}

	return dev;
}

/* Create program from a file and compile it */
cl_program build_program(cl_context ctx, cl_device_id dev, const char* filename) {

	cl_program program;
	FILE *program_handle;
	char *program_buffer, *program_log;
	size_t program_size, log_size;
	int err;

	/* Read program file and place content into buffer */
	program_handle = fopen(filename, "r");
	if (program_handle == NULL) {
		perror("Couldn't find the program file");
		exit(1);
	}
	fseek(program_handle, 0, SEEK_END);
	program_size = ftell(program_handle);
	rewind(program_handle);
	program_buffer = (char*) malloc(program_size + 1);
	program_buffer[program_size] = '\0';
	fread(program_buffer, sizeof(char), program_size, program_handle);
	fclose(program_handle);

	/* Create program from file */
	program = clCreateProgramWithSource(ctx, 1, (const char**) &program_buffer,
			&program_size, &err);
	if (err < 0) {
		perror("Couldn't create the program");
		exit(1);
	}
	free(program_buffer);

	/* Build program */
	err = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
	if (err < 0) {

		/* Find size of log and print to std output */
		clGetProgramBuildInfo(program, dev, CL_PROGRAM_BUILD_LOG, 0, NULL,
				&log_size);
		program_log = (char*) malloc(log_size + 1);
		program_log[log_size] = '\0';
		clGetProgramBuildInfo(program, dev, CL_PROGRAM_BUILD_LOG, log_size + 1,
				program_log, NULL);
		printf("%s\n", program_log);
		free(program_log);
		exit(1);
	}

	return program;
}

void initialize_cl() {

	/* Error var */
	cl_int err;

	/* Create device and context */
	printf("Getting OpenCL device ID... \n");
	cl_device_id device = create_device();
	printf("Creating the OpenCL context... \n");
	context = clCreateContext(NULL, 1, &device, NULL, NULL, &err);
	if (err < 0) {
		perror("Couldn't create a context");
		exit(1);
	}

	/* Build programs */
	program = build_program(context, device, PROGRAM_FILE_PATH);

	/* Create a command queue */
	printf("Creating OpenCL command queue... \n");
	queue = clCreateCommandQueue(context, device, 0, &err);
	if (err < 0) {
		perror("Couldn't create a command queue");
		exit(1);
	};

	// Create the kernels
	printf("Creating the OpenCL kernel... \n");
	kernel = clCreateKernel(program, "regression_analysis", &err);
	if (err < 0) {
		perror("Couldn't create a kernel");
		exit(1);
	};
}

int main() {

	printf("Initializing... \n");
	initialize_cl();
	printf("Done!\n");

	int mi, mj, mk;
	int count = 0;
	int size_x_rdcd = 0;
	int size_xtxinvxt_rdcd = 0;
	int voxels = 4;
	int nt = 450;
	int ixyz_bot = 0;
	int ixyz_top = voxels;

	/* Variables */
	int N = 450;
	int p = 14;
	int q = 12;
	int num_stimts = 8;
	float rms_min = 0.0;
	float min_lag[8] = {0, 0, 0, 0, 0, 0, 0, 0};
	float max_lag[8] = {0, 0, 0, 0, 0, 0, 0, 0};
	double * y_elts_buf = (double *) malloc(sizeof(double) * N * voxels);
	for(mi = 0; mi<450*voxels; ++mi) y_elts_buf[mi] = 0.9;
	int x_full_rows = 450;
	int x_full_cols = 14;
	double * x_full_elts = (double *)malloc(sizeof(double)*450*14);
	for(mi = 0; mi<450*14; ++mi) x_full_elts[mi] = 0.8;
	int xtxinv_full_rows = 14;
	int xtxinv_full_cols = 14;
	double * xtxinv_full_elts = (double *)malloc(sizeof(double)*14*14);
	for(mi = 0; mi<14*14; ++mi) xtxinv_full_elts[mi] = 0.7;
	int xtxinvxt_full_rows = 14;
	int xtxinvxt_full_cols = 450;
	double * xtxinvxt_full_elts = (double *)malloc(sizeof(double)*14*450);
	for(mi = 0; mi<14*450; ++mi) xtxinvxt_full_elts[mi] = 0.6;
	int x_base_rows = 450;
	int x_base_cols = 12;
	double * x_base_elts = (double *)malloc(sizeof(double)*450*12);
	for(mi = 0; mi<450*12; ++mi) x_base_elts[mi] = 0.5;
	int xtxinvxt_base_rows = 12;
	int xtxinvxt_base_cols = 450;
	double * xtxinvxt_base_elts = (double *)malloc(sizeof(double)*12*450);
	for(mi = 0; mi<12*450; ++mi) xtxinvxt_base_elts[mi] = (double)0.4;
	int * dims_x_rdcd = (int *)malloc(sizeof(int)*2*num_stimts);
	int * dims_xtxinvxt_rdcd = (int *)malloc(sizeof(int)*2*num_stimts);
	for(mi = 0; mi < num_stimts; ++mi){
		dims_x_rdcd[mi+count] = 450;
		dims_x_rdcd[mi+count+1] = 13;
		size_x_rdcd += 450*13;
		dims_xtxinvxt_rdcd[mi+count] = 13;
		dims_xtxinvxt_rdcd[mi+count+1] = 450;
		size_xtxinvxt_rdcd += 13*450;
		++count;
	}
	double * x_rdcd_1d = (double *)malloc(sizeof(double)*size_x_rdcd);
	for(mi = 0; mi<size_x_rdcd; ++mi) x_rdcd_1d[mi] = 0.3;
	double * xtxinvxt_rdcd_1d = (double *)malloc(sizeof(double)*size_xtxinvxt_rdcd);
	for(mi = 0; mi<size_xtxinvxt_rdcd; ++mi) xtxinvxt_rdcd_1d[mi] = 0.2;

	cl_int err = 0; //Error variable

	err = clSetKernelArg(kernel, 0, sizeof(int), &N); //Set the N arg
	err |= clSetKernelArg(kernel, 1, sizeof(int), &p); //Set the p arg
	err |= clSetKernelArg(kernel, 2, sizeof(int), &q); //Set the q arg
	err |= clSetKernelArg(kernel, 3, sizeof(int), &num_stimts); //Set up num_stimts arg
	if (err < 0) {perror("Couldn't set the N, p, q, and num_stimts kernel args!"); exit(1);}

	d_mat_1 = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double) * x_full_rows * x_full_cols, x_full_elts, &err); //Set up x_full.elts buf
	d_mat_2 = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double) * xtxinv_full_rows * xtxinv_full_cols, xtxinv_full_elts, &err); //Set up xtxinv_full.elts buf
	d_mat_3 = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double) * xtxinvxt_full_rows * xtxinvxt_full_cols, xtxinvxt_full_elts, &err); //Set up xtxinvxt_full.elts buf
	d_mat_4 = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double) * x_base_rows * x_base_cols, x_base_elts, &err); //Set up x_base.elts buf
	d_mat_5 = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double) * xtxinvxt_base_rows * xtxinvxt_base_cols, xtxinvxt_base_elts, &err); //Set up xtxinvxt_base.elts buf
	if (err < 0) {perror("Couldn't create and write to the x_full, xtxinv_full, xtxinvxt_full, x_base, and xtxinvxt_base cl buffers!"); exit(1);}
	err = clSetKernelArg(kernel, 4, sizeof(int), &x_full_rows); //Set the x_full.rows arg
	err |= clSetKernelArg(kernel, 5, sizeof(int), &x_full_cols); //Set the x_full.cols arg
	err |= clSetKernelArg(kernel, 6, sizeof(cl_mem), &d_mat_1); //Set the x_full.elts arg
	err |= clSetKernelArg(kernel, 7, sizeof(int), &xtxinv_full_rows); //Set the xtxinv_full.rows arg
	err |= clSetKernelArg(kernel, 8, sizeof(int), &xtxinv_full_cols); //Set the xtxinv_full.cols arg
	err |= clSetKernelArg(kernel, 9, sizeof(cl_mem), &d_mat_2); //Set the xtxinv_full.elts arg
	err |= clSetKernelArg(kernel, 10, sizeof(int), &xtxinvxt_full_rows); //Set the xtxinvxt_full.rows arg
	err |= clSetKernelArg(kernel, 11, sizeof(int), &xtxinvxt_full_cols); //Set the xtxinvxt_full.cols arg
	err |= clSetKernelArg(kernel, 12, sizeof(cl_mem), &d_mat_3); //Set the xtxinvxt_full.elts arg
	err |= clSetKernelArg(kernel, 13, sizeof(int), &x_base_rows); //Set the x_base.rows arg
	err |= clSetKernelArg(kernel, 14, sizeof(int), &x_base_cols); //Set the x_base.cols arg
	err |= clSetKernelArg(kernel, 15, sizeof(cl_mem), &d_mat_4); //Set the x_base.elts arg
	err |= clSetKernelArg(kernel, 16, sizeof(int), &xtxinvxt_base_rows); //Set the xtxinvxt_base.rows arg
	err |= clSetKernelArg(kernel, 17, sizeof(int), &xtxinvxt_base_cols); //Set the xtxinvxt_base.cols arg
	err |= clSetKernelArg(kernel, 18, sizeof(cl_mem), &d_mat_5); //Set the xtxinvxt_base.elts arg
	if (err < 0) {perror("Couldn't set the x_full, xtxinv_full, xtxinvxt_full, x_base, and xtxinvxt_base kernel args!"); exit(1);}

	int_buf_1 = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(int) * 2 * num_stimts, dims_x_rdcd, &err); //for i x_rdcd
	d_mat_ptr_1 = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double) * size_x_rdcd, x_rdcd_1d, &err); //for i x_rdcd
	int_buf_2 = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(int) * 2 * num_stimts, dims_xtxinvxt_rdcd, &err); //for i xtxinvxt_rdcd
	d_mat_ptr_2 = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double) * size_xtxinvxt_rdcd, xtxinvxt_rdcd_1d, &err); //for i xtxinvxt_rdcd
	if (err < 0) {perror("Couldn't create and write to the x_rdcd and xtxinvxt_rdcd cl buffers!"); exit(1);}
	err = clSetKernelArg(kernel, 19, sizeof(cl_mem), &int_buf_1); //Set the combined x_rdcd[].rows and x_rdcd[].cols arg
	err |= clSetKernelArg(kernel, 20, sizeof(cl_mem), &d_mat_ptr_1); //Set the unpacked x_rdcd[].elts arg
	err |= clSetKernelArg(kernel, 21, sizeof(cl_mem), &int_buf_2); //Set the combined xtxinvxt_rdcd[].rows and xtxinvxt_rdcd[].cols arg
	err |= clSetKernelArg(kernel, 22, sizeof(cl_mem), &d_mat_ptr_2); //Set the unpacked xtxinvxt_rdcd[].elts arg
	if (err < 0) {perror("Couldn't set the x_rdcd and xtxinvxt_rdcd kernel args!"); exit(1);}

	d_vec_1 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(double) * N * voxels, y_elts_buf, &err); //Set up y_elts_buf buf
	if (err < 0) {perror("Couldn't create and write to the y_elts_buf cl buffer!"); exit(1);}
	err = clSetKernelArg(kernel, 23, sizeof(cl_mem), &d_vec_1);
	if (err < 0) {perror("Couldn't set the y_elts_buf kernel arg!"); exit(1);}

	err = clSetKernelArg(kernel, 24, sizeof(float), &rms_min); //Set the rms_min arg
	if (err < 0) {perror("Couldn't set the rms_min kernel arg!"); exit(1);}
	//END INPUTS

	//OUTPUTS
	//Need to malloc memory for all of them
	float * mse_buf = (float *)calloc(voxels, sizeof(float));
	double * coef_elts = (double *)calloc(p * voxels, sizeof(double));
	double * scoef_elts = (double *)calloc(p * voxels, sizeof(double));
	double * tcoef_elts = (double *)calloc(p * voxels, sizeof(double));
	float * fpart_buffer = (float *)calloc(num_stimts * voxels, sizeof(float));
	float * rpart_buffer = (float *)calloc(num_stimts * voxels, sizeof(float));
	float * ffull_buf= (float *)calloc(voxels, sizeof(float));
	float  * rfull_buf= (float *)calloc(voxels, sizeof(float));
	int  * novar_buf = (int *)calloc(voxels, sizeof(int));
	float * fitts_buf = (float *)calloc(nt * voxels, sizeof(float));
	float * errts_buf = (float *)calloc(nt * voxels, sizeof(float));

	return_int_ptr_1 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(float) * voxels, mse_buf, &err); //Set up mse output buf
	if (err < 0) {perror("Couldn't create the mse cl output buffer!"); exit(1);}
	err = clSetKernelArg(kernel, 25, sizeof(cl_mem), &return_int_ptr_1); //Set the mse arg
	if (err < 0) {perror("Couldn't set the mse kernel arg!"); exit(1);}

	return_vec_1 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(double) * p * voxels, coef_elts, &err); //Set up coef.elts buf
	return_vec_2 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(double) * p * voxels, scoef_elts, &err); //Set up scoef.elts buf
	return_vec_3 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(double) * p * voxels, tcoef_elts, &err); //Set up tcoef.elts buf
	if (err < 0) {perror("Couldn't create the coef.elts, scoef.elts, and tcoef.elts cl output buffers!"); exit(1);}
	err = clSetKernelArg(kernel, 26, sizeof(cl_mem), &return_vec_1); //Set the coef.elts arg
	err |= clSetKernelArg(kernel, 27, sizeof(cl_mem), &return_vec_2); //Set the scoef.elts arg
	err |= clSetKernelArg(kernel, 28, sizeof(cl_mem), &return_vec_3); //Set the tcoef.elts arg
	if (err < 0) {printf("Couldn't set the coef, scoef, and tcoef kernel args! %i", err); exit(1);}

	return_flt_buf_1 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(float) * num_stimts * voxels, fpart_buffer, &err); //for i/o fpart
	return_flt_buf_2 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(float) * num_stimts * voxels, rpart_buffer, &err); //for i/o rpart
	if (err < 0) {printf("Couldn't create the fpart and rpart cl buffers! %i", err); exit(1);}
	err = clSetKernelArg(kernel, 29, sizeof(cl_mem), &return_flt_buf_1);
	err |= clSetKernelArg(kernel, 30, sizeof(cl_mem), &return_flt_buf_2);
	if (err < 0) {perror("Couldn't set the fpart and rpart kernel args!"); exit(1);}

	return_flt_ptr_1 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(float) * voxels, ffull_buf, &err); //Set up ffull output buf
	return_flt_ptr_2 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(float) * voxels, rfull_buf, &err); //Set up rfull output buf
	return_int_ptr_2 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(int) * voxels, novar_buf, &err); //Set up novar output buf
	if (err < 0) {perror("Couldn't create the ffull, rfull, and novar cl output buffers!"); exit(1);}
	err = clSetKernelArg(kernel, 31, sizeof(cl_mem), &return_flt_ptr_1);
	err |= clSetKernelArg(kernel, 32, sizeof(cl_mem), &return_flt_ptr_2);
	err |= clSetKernelArg(kernel, 33, sizeof(cl_mem), &return_int_ptr_2);
	if (err < 0) {perror("Couldn't set the ffull, rfull, and novar kernel args!"); exit(1);}

	return_flt_buf_3 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(float) * nt * voxels, fitts_buf, &err); //Set up fitts output buf
	return_flt_buf_4 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(float) * nt * voxels, errts_buf, &err); //Set up errts output buf
	if (err < 0) {perror("Couldn't create the errts and fitts cl buffers!"); exit(1);}
	err = clSetKernelArg(kernel, 34, sizeof(cl_mem), &return_flt_buf_3);
	err |= clSetKernelArg(kernel, 35, sizeof(cl_mem), &return_flt_buf_4);
	if (err < 0) {perror("Couldn't set the fitts_buf and errts_buf kernel args!"); exit(1);}
	//END OUTPUTS

	//Oops, I'm lazy
	//Malloc memory for the temp bufs
	double * coef_temp_buf = (double *)calloc(p * voxels, sizeof(double) );
	double * yhat_buf = (double *)calloc(N * voxels, sizeof(double) );
	double * e_buf = (double *)calloc(N * voxels, sizeof(double) );

	int_buf_3 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(int) * num_stimts, min_lag, &err); //Set up min_lag buf
	int_buf_4 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(int) * num_stimts, max_lag, &err); //Set up max_lag buf
	if (err < 0) {perror("Couldn't create the min_lag and max_lag cl buffers!"); exit(1);}
	err = clSetKernelArg(kernel, 36, sizeof(cl_mem), &int_buf_3);
	err |= clSetKernelArg(kernel, 37, sizeof(cl_mem), &int_buf_4);
	if (err < 0) {perror("Couldn't set the min_lag and max_lag kernel args!"); exit(1);}
	d_vec_1 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(double) * p * voxels, coef_temp_buf, &err); //Set up coef_temp_buf buf
	if (err < 0) {perror("Couldn't create the coef_temp_buf cl buffer!"); exit(1);}
	err = clSetKernelArg(kernel, 38, sizeof(cl_mem), &d_vec_1);
	if (err < 0) {perror("Couldn't set the coef_temp_buf arg!"); exit(1);}
	d_vec_2 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(double) * N * voxels, yhat_buf, &err); //Set up yhat_buf buf
	d_vec_3 = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(double) * N * voxels, e_buf, &err); //Set up e_buf buf
	if (err < 0) {perror("Couldn't create the yhat_buf and e_buf cl buffer!"); exit(1);}
	err = clSetKernelArg(kernel, 39, sizeof(cl_mem), &d_vec_2);
	err |= clSetKernelArg(kernel, 40, sizeof(cl_mem), &d_vec_3);
	if (err < 0) {perror("Couldn't set the yhat_buf and e_buf arg!"); exit(1);}
	//END LAZY

	err = clSetKernelArg(kernel, 41, sizeof(int), &ixyz_bot);
	err |= clSetKernelArg(kernel, 42, sizeof(int), &ixyz_top);
	if (err < 0) {perror("Couldn't set the ixyz_bot and ixyz_top kernel arg!"); exit(1);}

	printf("We did it, all arguments were set successfully! Enqueue the kernels...\n");

	//Set the global workgroup size; 32 is used because the Intel Phi likes it
	//No local workgroup size is define; OpenCL will pick what is best
	size_t global_size = voxels;

	/*int * temp_buf = (int *)calloc(voxels, sizeof(int));

	return_temp_buf = clCreateBuffer(context,
					CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
					sizeof(int) * voxels, temp_buf, &err); //Set up errts output buf
	if (err < 0) {
		perror("Couldn't create the temp cl buffer!");
		exit(1);
	}

	err = clSetKernelArg(kernel, 43, sizeof(cl_mem), &return_temp_buf);
	if (err < 0) {
		printf("Couldn't set the temp kernel arg! %i", err);
		exit(1);
	}*/

	cl_event event = NULL;

	/* Enqueue kernel */
	err = clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &global_size, NULL, 0, NULL, &event);
	if (err < 0) {
		printf("Couldn't enqueue the kernel, err num %i\n", err);
		exit(1);
	}
	printf("Kernels enqueued successfully! Waiting for them to finish...\n");

	err = clWaitForEvents(1, &event);

	int status = clReleaseEvent(event);

	/*err = clEnqueueReadBuffer(queue, return_temp_buf, CL_TRUE, 0,
					sizeof(int) * voxels, temp_buf, 0, NULL, NULL);
	if (err < 0) {
		perror("Couldn't read the temp buffer!");
		exit(1);
	}

	for(mi = 0; mi < voxels; ++mi) printf("%i ,", temp_buf[mi]);
	printf("\n");*/

	/* Read those outputs */
	printf("Reading mse_buf...\nmse: ");
	err = clEnqueueReadBuffer(queue, return_int_ptr_1, CL_TRUE, 0, sizeof(float) * voxels, mse_buf, 0, NULL, NULL);
	if (err < 0) {printf("Couldn't read the mse cl buffer! %i", err); exit(1);}
	for(mi = 0; mi<voxels; ++mi) printf("%f, ", mse_buf[mi]);

	printf("\nReading coef_elts...\ncoef: ");
	err = clEnqueueReadBuffer(queue, return_vec_1, CL_TRUE, 0, sizeof(double) * p * voxels, coef_elts, 0, NULL, NULL);
	for(mi = 0; mi<p*voxels; ++mi) printf("%f, ", coef_elts[mi]);
	printf("\nReading scoef_elts...\nscoef: ");
	err |= clEnqueueReadBuffer(queue, return_vec_1, CL_TRUE, 0, sizeof(double) * p * voxels, scoef_elts, 0, NULL, NULL);
	for(mi = 0; mi<p*voxels; ++mi) printf("%f, ", scoef_elts[mi]);
	printf("\nReading tcoef_elts...\ntcoef: ");
	err |= clEnqueueReadBuffer(queue, return_vec_1, CL_TRUE, 0, sizeof(double) * p * voxels, tcoef_elts, 0, NULL, NULL);
	for(mi = 0; mi<p*voxels; ++mi) printf("%f, ", tcoef_elts[mi]);
	if (err < 0) {perror("Couldn't read the coef, scoef, and tcoef cl buffers!"); exit(1);}

	printf("\nReading fpart_buf...\nfpart: ");
	err = clEnqueueReadBuffer(queue, return_flt_buf_1, CL_TRUE, 0, sizeof(float) * num_stimts * voxels, fpart_buffer, 0, NULL, NULL);
	for(mi = 0; mi<num_stimts*voxels; ++mi) printf("%f, ", fpart_buffer[mi]);
	printf("\nReading rpart_buf...\nrpart: ");
	err |= clEnqueueReadBuffer(queue, return_flt_buf_2, CL_TRUE, 0, sizeof(float) * num_stimts * voxels, rpart_buffer, 0, NULL, NULL);
	for(mi = 0; mi<num_stimts*voxels; ++mi) printf("%f, ", rpart_buffer[mi]);
	if (err < 0) {perror("Couldn't read the fpart and rpart cl buffers!"); exit(1);}

	printf("\nReading ffull_buf...\nffull: ");
	err = clEnqueueReadBuffer(queue, return_flt_ptr_1, CL_TRUE, 0, sizeof(float) * voxels, ffull_buf, 0, NULL, NULL);
	for(mi = 0; mi<voxels; ++mi) printf("%f, ", ffull_buf[mi]);
	printf("\nReading rfull_buf...\nrfull: ");
	err |= clEnqueueReadBuffer(queue, return_flt_ptr_2, CL_TRUE, 0, sizeof(float) * voxels, rfull_buf, 0, NULL, NULL);
	for(mi = 0; mi<voxels; ++mi) printf("%f, ", rfull_buf[mi]);
	printf("\nReading novar_buf...\nnovar: ");
	err |= clEnqueueReadBuffer(queue, return_int_ptr_2, CL_TRUE, 0, sizeof(int) * voxels, novar_buf, 0, NULL, NULL);
	for(mi = 0; mi<voxels; ++mi) printf("%f, ", novar_buf[mi]);
	if (err < 0) {perror("Couldn't read the ffull, rfull, and novar cl buffers!"); exit(1);}

	printf("\nReading fitts...\nfitts: ");
	err = clEnqueueReadBuffer(queue, return_flt_buf_3, CL_TRUE, 0, sizeof(float) * nt * voxels, fitts_buf, 0, NULL, NULL);
	for(mi = 0; mi<nt*voxels; ++mi) printf("%f, ", fitts_buf[mi]);
	printf("\nReading errts...\nerrts: ");
	err |= clEnqueueReadBuffer(queue, return_flt_buf_4, CL_TRUE, 0, sizeof(float) * nt * voxels, errts_buf, 0, NULL, NULL);
	for(mi = 0; mi<nt*voxels; ++mi) printf("%f, ", errts_buf[mi]);
	if (err < 0) {perror("Couldn't read the fitts_buf and errts_buf cl buffers!"); exit(1);}

	printf("\nAll kernels finished! Releasing input memory...\n");

	  /* Deallocate OpenCL input resources */
	  free (y_elts_buf);
	  free (dims_x_rdcd);
	  free (dims_xtxinvxt_rdcd);
	  free (x_rdcd_1d);
	  free (xtxinvxt_rdcd_1d);
	  free(x_full_elts);
	  free(xtxinv_full_elts);
	  free(xtxinvxt_full_elts);
	  free(x_base_elts);
	  free(xtxinvxt_base_elts);
	  free(coef_temp_buf);
	  free(yhat_buf);
	  free(e_buf);
	  clReleaseCommandQueue(queue);
	  clReleaseProgram(program);
	  clReleaseContext(context);
	  clReleaseMemObject(int_buf_1);
	  clReleaseMemObject(int_buf_2);
	  clReleaseMemObject(int_buf_3);
	  clReleaseMemObject(int_buf_4);
	  clReleaseMemObject(d_mat_1);
	  clReleaseMemObject(d_mat_2);
	  clReleaseMemObject(d_mat_3);
	  clReleaseMemObject(d_mat_4);
	  clReleaseMemObject(d_mat_5);
	  clReleaseMemObject(d_mat_ptr_1);
	  clReleaseMemObject(d_mat_ptr_2);
	  clReleaseMemObject(d_vec_1);
	  clReleaseMemObject(return_flt_buf_1);
	  clReleaseMemObject(return_flt_buf_2);
	  clReleaseMemObject(return_flt_buf_3);
	  clReleaseMemObject(return_flt_buf_4);
	  clReleaseMemObject(return_int_ptr_1);
	  clReleaseMemObject(return_int_ptr_2);
	  clReleaseMemObject(return_vec_1);
	  clReleaseMemObject(return_vec_2);
	  clReleaseMemObject(return_vec_3);
	  clReleaseMemObject(return_flt_ptr_1);
	  clReleaseMemObject(return_flt_ptr_2);
	  clReleaseKernel(kernel);
		free(mse_buf);
		free(coef_elts);
		free(scoef_elts);
		free(tcoef_elts);
		free(fpart_buffer);
		free(rpart_buffer);
		free(ffull_buf);
		free(rfull_buf);
		free(novar_buf);
		free(fitts_buf);
		free(errts_buf);

	printf("Exiting...");

	return (0);
}