 	CUDA_SAFE_CALL (cudaMemcpy (b_d, b, size * sizeof(Real), cudaMemcpyHostToDevice));
 	CUDA_SAFE_CALL (cudaMemcpy (temp_d, temp, size_temp * sizeof(Real), cudaMemcpyHostToDevice));
 	
-	// block size
-	const uint block_size = 64;
-	
 	// block and grid dimensions
 	
 	///////////////////////////////////////

	dim3 dimGrid (num_rows - 2, (num_cols - 2) / block_size);
	///////////////////////////////////////
	
		
	// iteration loop
	for (iter = 1; iter <= it_max; ++iter) {
		
		// update red cells
		red_kernel <<<dimGrid, dimBlock>>> (num_rows, num_cols, aP_d, aW_d, aE_d, aS_d, aN_d, b_d, temp_d);
		
		// update black cells
		black_kernel <<<dimGrid, dimBlock>>> (num_rows, num_cols, aP_d, aW_d, aE_d, aS_d, aN_d, b_d, temp_d);
		
		// sync threads (needed?)
		CUDA_SAFE_CALL (cudaThreadSynchronize());
		
		// transfer memory back to host to check for convergence
		CUDA_SAFE_CALL (cudaMemcpy (temp, temp_d, size_temp * sizeof(Real), cudaMemcpyDeviceToHost));
		
		// check residual to see if done with iterations
		Real norm_L2 = 0.0;
		for (uint i = 0; i < size_temp; ++i) {
			norm_L2 += (temp[i] - temp_old[i]) * (temp[i] - temp_old[i]);
			temp_old[i] = temp[i];
		}
		norm_L2 = sqrt(norm_L2 / (size));
		
		// if tolerance has been reached, end SOR iterations
		if (norm_L2 < tol) {
			break;
		}	
	}
	
	// free device memory
	CUDA_SAFE_CALL (cudaFree(aP_d));

										   + aS[ind] * temp[col * rowmax + (row - 1)]
										   + aN[ind] * temp[col * rowmax + (row + 1)]);
	
		//temp[ind2] = temp[ind2] * (1.0 - omega) + omega * (res / aP[ind]);
		Real temp_old = temp[ind2];
		Real temp_new = temp_old * (1.0 - omega) + omega * (res / aP[ind]);
		
		temp[ind2] = temp_new;
		res = temp_old - temp_new;
		
		// store squared residual from each thread
		res_cache[threadIdx.y] = res * res;
		
		// synchronize threads in block
		__syncthreads();
		
		// add up squared residuals for block
		uint i = block_size / 2;
		while (i != 0) {
			if (threadIdx.y < i) {
				res_cache[threadIdx.y] += res_cache[threadIdx.y + i];
			}
			__syncthreads();
			i /= 2;
		}
		
		// store block's summed residuals
		if (threadIdx.y == 0) {
			bl_norm_L2[blockIdx.x + (gridDim.x * blockIdx.y)] = res_cache[0];
		}
	}
}

__global__ void black_kernel (uint rowmax, uint colmax, const Real * aP,
														  const Real * aW, const Real * aE, const Real * aS,
														  const Real * aN, const Real * b, Real * temp, Real * bl_norm_L2)
{	
	uint row = 1 + (blockIdx.y * blockDim.y) + threadIdx.y;
	uint col = 1 + (blockIdx.x * blockDim.x) + threadIdx.x;
	
	// store residual for block
	__shared__ Real res_cache[block_size];
	
	res_cache[threadIdx.y] = 0.0;
	
	// only black cell if odd
	if ((row + col) % 2 == 1) {
