 		// Calculate the data costs on the other levels
 		for (int t = 1; t < numLevels; t++) {
 			// Calculate the dimensions of the grid at this level
-			int w = (int)ceil(data[t - 1]->getSize(0) / 2.0), h = (int)ceil(data[t - 1]->getSize(1) / 2.0);
+			int w = widths[t] = (int)ceil(widths[t - 1] / 2.0);
+			int h = heights[t] = (int)ceil(heights[t - 1] / 2.0);
 
 			// We gotta make sure that we don't get too 'coarse'
-			if (w < 1 || h < 1) return;
-
-			// The data cost of each node GINF_BP_BLOCK_WIDTHon the current level uses the sum of
-			// the data costs of four nodes below this level
-			data[t] = new Matrix<T>(3, w, h, grid->getNumLabels());
-			for (int y = 0; y < data[t - 1]->getSize(1); y++) {
-				for (int x = 0; x < data[t - 1]->getSize(0); x++) {
-					for (int i = 0; i < grid->getNumLabels(); i++) {
-						data[t]->get(x / 2, y / 2, i) += data[t - 1]->get(x, y, i);
-					}
-				}
-			}
-			
-			// Copy data cost matrix on the device
-			dData[t] = gpu.createRawMat(data[t]);
-			//gpuBpCalcDataCost<<<gridSize, blockSize>>>(dGrid, w, h, data[t - 1]->getSize(0), data[t - 1]->getSize(1), dData[t], dData[t - 1]);
+			if (w < 1 || h < 1) {
+				numLevels = t - 1;
+				break;
+			}
+
+			// The data cost of each node on the current level uses the sum of the 
+			// data costs of four nodes below this level. We calculate this in a 
+			dData[t] = gpu.createRawMat<T>(w * h * grid->getNumLabels());
+			gpuBpCalcDataCost<<<gridSize, blockSize>>>(dGrid, w, h, widths[t - 1], heights[t - 1], dData[t], dData[t - 1]);
 		}
 		
 		// For each level, we'll do message passing
 		for (int t = numLevels - 1; t >= 0; t--) {
-			int w = data[t]->getSize(0), h = data[t]->getSize(1);
+			int w = widths[t], h = heights[t];
 			
 			// Allocate message matrices on the device for this particular level
 			dMsgs[t] = gpu.createRawMat<T>(w * h * GINF_NUM_DIR * grid->getNumLabels());
 			dFrom[t] = gpu.createRawMat<T>(w * h * grid->getNumLabels());
+			
+			// Calculate gridSizes
+			dim3 lvlGridSize(GINF_DIM(w, blockSize.x), GINF_DIM(h, blockSize.y));
+			dim3 lvlGridSizeBipartite(GINF_DIM(w / 2 + 1, blockSize.x), GINF_DIM(h, blockSize.y));
 
 			// If we're not on the top level, the messages are initialized from the previous level
 			if (t != numLevels - 1) {
-				int pw = data[t + 1]->getSize(0), ph = data[t + 1]->getSize(1);
-				
-				gpuBpInitMsgs<<<gridSize, blockSize>>>(dGrid, w, h, pw, ph, dMsgs[t], dMsgs[t + 1]);
-
-				// Free up memory
-				delete data[t + 1];
+				gpuBpInitMsgs<<<lvlGridSize, blockSize>>>(dGrid, w, h, widths[t + 1], heights[t + 1], dMsgs[t], dMsgs[t + 1]);
 			}
 			
 			// For each iteration, run belief propagation
-			dim3 lvlGridSize(GINF_DIM(w / 2 + 1, blockSize.x), GINF_DIM(h, blockSize.y));
 			for (int i = 0; i < numItersPerLevel; i++) {
-				gpuBpCollectKernel<<<gridSize, blockSize>>>(dGrid, w, h, dMsgs[t], dFrom[t], dData[t]);
+				gpuBpCollectKernel<<<lvlGridSize, blockSize>>>(dGrid, w, h, dMsgs[t], dFrom[t], dData[t]);
 
 				// For each direction, send a message
 				for (int d = 0; d < GINF_NUM_DIR; d++) {
-					gpuBpKernel<<<lvlGridSize, blockSize>>>(i, dGrid, w, h, dMsgs[t], dFrom[t], dData[t], d);
+					gpuBpKernel<<<lvlGridSizeBipartite, blockSize>>>(i, dGrid, w, h, dMsgs[t], dFrom[t], dData[t], d);
 				}
 			}
 		}

		}
	}
	
	/*
	template <typename T>
	__global__ void gpuBpCalcDataCost(GpuGrid<T> *grid, int w, int h, int pw, int ph, T *data, T *prev) {
		GINF_DECL_X; GINF_DECL_Y;
		MatDim dim(w, h, grid->getNumLabels());
		MatDim dimPrev(pw, ph, grid->getNumLabels());
		
		if (dim.isValid(x / 2, y / 2)) {
			for (int i = 0; i < grid->getNumLabels(); i++) {
				data[dim.idx(x / 2, y / 2, i)] += prev[dimPrev.idx(x, y, i)];
			}
		}
	}*/
	
	template <typename T>
	__global__ void gpuBpInitMsgs(GpuGrid<T> *grid, int w, int h, int pw, int ph, T *msgs, T *prev) {
		GINF_DECL_X; GINF_DECL_Y;	
		MatDim dim(w, h, GINF_NUM_DIR);
		MatDim dimPrev(pw, ph, GINF_NUM_DIR);
		
		if (dim.isValidStrict(x, y)) {
			for (int d = 0; d < GINF_NUM_DIR; d++) {
				for (int i = 0; i < grid->getNumLabels(); i++) {
					// Initialize message with node on the above level
					msgs[dim(x, y, d, i)] = prev[dimPrev(x / 2, y / 2, d, i)];
				}
			}
		}

		GpuGrid<T> *dGrid = gpu.createGrid(grid);
		
		// At the finest level, the data costs correspond to the original problem
		widths[0] = grid->getWidth();
		heights[0] = grid->getHeight();

		// Copy the original data costs to the device
		dData[0] = gpu.createRawMat(grid->getDataCosts());

		// Prepare block sizes
		dim3 blockSize(GINF_BP_BLOCK_WIDTH, GINF_BP_BLOCK_HEIGHT);
