     __syncthreads();
 
     if (pattern < totalPatterns) { // Remove padded threads!
-        for(i = 0; i < PADDED_STATE_COUNT; i++) {
-            sum1 += sMatrix1[i * 4 + state] * sPartials1[patIdx * 16 + pat * 4 + i];
-            sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
-        }
+//        for(i = 0; i < PADDED_STATE_COUNT; i++) {
+//            sum1 += sMatrix1[i * 4 + state] * sPartials1[patIdx * 16 + pat * 4 + i];
+//            sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
+//        }
+
+//        sum1 = sMatrix1[state] * sPartials1[patIdx * 16 + pat * 4];
+//        sum2 = sMatrix2[state] * sPartials2[patIdx * 16 + pat * 4];
+//
+//        sum1 += sMatrix1[4 + state] * sPartials1[patIdx * 16 + pat * 4 + 1];
+//        sum2 += sMatrix2[4 + state] * sPartials2[patIdx * 16 + pat * 4 + 1];
+//
+//        sum1 += sMatrix1[8 + state] * sPartials1[patIdx * 16 + pat * 4 + 2];
+//        sum2 += sMatrix2[8 + state] * sPartials2[patIdx * 16 + pat * 4 + 2];
+//
+//        sum1 += sMatrix1[12 + state] * sPartials1[patIdx * 16 + pat * 4 + 3];
+//        sum2 += sMatrix2[12 + state] * sPartials2[patIdx * 16 + pat * 4 + 3];
+
+        i = state;
+        sum1 = sMatrix1[i * 4 + state] * sPartials1[patIdx * 16 + pat * 4 + i];
+        sum2 = sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
+
+        i = (++i) & 0x3;
+        sum1 += sMatrix1[i * 4 + state] * sPartials1[patIdx * 16 + pat * 4 + i];
+        sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
+
+        i = (++i) & 0x3;
+        sum1 += sMatrix1[i * 4 + state] * sPartials1[patIdx * 16 + pat * 4 + i];
+        sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
+
+        i = (++i) & 0x3;
+        sum1 += sMatrix1[i * 4 + state] * sPartials1[patIdx * 16 + pat * 4 + i];
+        sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
+
         partials3[u] = sum1 * sum2 / fixedScalingFactors[patIdx * 4 + pat];
     }
 

    int pattern = __umul24(blockIdx.x, PATTERN_BLOCK_SIZE * 4) + patIdx * 4 + pat;

//  int deltaPartialsByState = __umul24(pattern,PADDED_STATE_COUNT);
    int deltaPartialsByState = 4 * 4 * (blockIdx.x * PATTERN_BLOCK_SIZE + patIdx);
    int deltaPartialsByMatrix = __umul24(matrix, __umul24( PADDED_STATE_COUNT, patternCount));

    int x2 = __umul24(matrix, PADDED_STATE_COUNT * PADDED_STATE_COUNT);

    REAL* matrix1 = matrices1 + x2; // Points to *this* matrix
    REAL* matrix2 = matrices2 + x2;

                                                              int totalPatterns) {

    int tx = threadIdx.x;
    int state = tx & 0x3; //  tx % 4;
    int pat = tx / 4;
    int patIdx = threadIdx.y;
    int matrix = blockIdx.y;
