 		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) {
+//__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 kernelStatesStatesByPatternBlockCoherent(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;

	__shared__ REAL sMatrix1[16];
	__shared__ REAL sMatrix2[16];

//	__shared__ REAL sPartials1[PATTERN_BLOCK_SIZE*4*4];

	__shared__ INT sStates1[PATTERN_BLOCK_SIZE*4];
	__shared__ REAL sPartials2[PATTERN_BLOCK_SIZE*4*4];

	// copy PADDED_STATE_COUNT*PATTERN_BLOCK_SIZE lengthed partials
	if (pattern < totalPatterns) {
//		sPartials1[patIdx*16 + tx] = partials1[y + tx]; // All coalesced memory reads
		sPartials2[patIdx*16 + tx] = partials2[y + tx];
		if (patIdx < PATTERN_BLOCK_SIZE/4)
			sStates1[patIdx*16 + tx] = states1[blockIdx.x*PATTERN_BLOCK_SIZE/4 + tx];
	} else {
//		sPartials1[patIdx*16 + tx] = 0;
		sPartials2[patIdx*16 + tx] = 0;
	}

	if (patIdx == 0 ) {


}


__global__ void kernelStatesStatesByPatternBlockCoherentSmall(int* states1, int* states2,
		REAL* partials3, REAL* matrices1, REAL* matrices2, int totalPatterns) {

	int tx = threadIdx.x;
	int state = tx % 4;
