 	REAL *matrix1 = matrices1 + x2; // Points to *this* matrix
 	REAL *matrix2 = matrices2 + x2;
 
-
-	int y = deltaPartialsByState + deltaPartialsByMatrix;
 	int u = tx + deltaPartialsByState + deltaPartialsByMatrix;
 
 #ifdef KERNEL_PRINT_ENABLED

		partials3[u] = sum1 * sum2;
}

__global__ void kernelStatesPartials(INT* states1, REAL* partials2, REAL* partials3, REAL* matrices1, REAL* matrices2) {

	REAL sum = 0;

	int i;

	// blockIdx.y = matrix; blockIdx.x = pattern; threadIdx.x = state
	// blockDim.x = PADDED_STATE_COUNT; gridDim.y = matrixCount; gridDim.x = patternCount;

	int deltaPartialsByState = blockIdx.x * blockDim.x;
	int deltaPartialsByMatrix = blockIdx.y * blockDim.x * gridDim.x;

//	int x = threadIdx.x * blockDim.x + blockIdx.y * blockDim.x * blockDim.x;
	int y = deltaPartialsByState + deltaPartialsByMatrix;

	int matrixOffset = blockIdx.y * PADDED_STATE_COUNT * PADDED_STATE_COUNT +  threadIdx.x * PADDED_STATE_COUNT;

	int u = threadIdx.x + deltaPartialsByState + deltaPartialsByMatrix;

	int state1 = (int) states1[blockIdx.x]; // TODO Should cache in shared memory

#ifdef KERNEL_PRINT_ENABLED_SP
	printf("thread #(%d,%d,%d): x = %d, y = %d, u = %d, state = %d, ",blockIdx.x,blockIdx.y,threadIdx.x,x,y,u,state1); // for debugging in emulation_mode
#endif


	if (state1 < PADDED_STATE_COUNT) { // < PADDED_STATE_COUNT

		REAL value = matrices1[matrixOffset + state1*blockDim.x + threadIdx.x]; // TODO should cache value

		for(i=0; i<PADDED_STATE_COUNT; i++) {
			sum += matrices2[matrixOffset + i] * partials2[y + i]; // TODO should block read
		}

#ifdef KERNEL_PRINT_ENABLED_SP
		printf("sum = %1.2e, value = %1.2e\n",sum,value);
#endif

		partials3[u] = sum * value;

	} else { // Missing data

		for(i=0; i<PADDED_STATE_COUNT; i++) {
			sum += matrices2[matrixOffset + i] * partials2[y + i];
		}

		partials3[u] = sum;
	}
}


__global__ void kernelStatesStates(INT* states1, INT* states2, REAL* partials3, REAL* matrices1, REAL* matrices2) {

	// blockIdx.y = matrix; blockIdx.x = pattern; threadIdx.x = state
	// blockDim.x = PADDED_STATE_COUNT; gridDim.y = matrixCount; gridDim.x = patternCount;

	INT* states1, INT* states2, REAL* partials3, REAL* matrices1, REAL* matrices2, REAL* scalingFactors,
	const unsigned int patternCount, const unsigned int matrixCount, int doRescaling) {

#if (PADDED_STATE_COUNT == 4)
	dim3 block(16,PATTERN_BLOCK_SIZE);
	dim3 grid(patternCount/(PATTERN_BLOCK_SIZE*4), matrixCount);
	if (patternCount % (PATTERN_BLOCK_SIZE*4) != 0)
		grid.x +=1;

	kernelStatesStatesByPatternBlockCoherentSmall<<<grid, block>>>(states1, states2, partials3, matrices1, matrices2, patternCount);
#else
//	dim3 grid(patternCount/PATTERN_BLOCK_SIZE, matrixCount);
//	if (patternCount % PATTERN_BLOCK_SIZE != 0)
//		grid.x += 1;
//	dim3 block(PADDED_STATE_COUNT,PATTERN_BLOCK_SIZE);
//
//	kernelStatesStatesByPatternBlockCoherent<<<grid, block>>>(states1, states2, partials3, matrices1, matrices2, patternCount);

	// Defunct
	dim3 grid(patternCount, matrixCount);
	dim3 block(PADDED_STATE_COUNT,1);

	kernelStatesStates<<<grid, block>>>(states1, states2, partials3, matrices1, matrices2);
#endif

#ifdef DEBUG
	fprintf(stderr,"Completed GPU SS\n");
