 
         if (state1 < PADDED_STATE_COUNT)
             sum1 = sMatrix1[state1 * 4 + state];
-        else
-            sum1 = 1.0;
-
-        for(i=0; i<PADDED_STATE_COUNT; i++) {
-            sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
-        }
+            // sum1 is initialized to 1
+//        else
+//            sum1 = 1.0;
+
+//        for(i=0; i<PADDED_STATE_COUNT; i++) {
+//            sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
+//        }
+
+        i = state;
+        sum2 = sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
+        i = (++i) & 0x3;
+        sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
+        i = (++i) & 0x3;
+        sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
+        i = (++i) & 0x3;
+        sum2 += sMatrix2[i * 4 + state] * sPartials2[patIdx * 16 + pat * 4 + i];
         partials3[u] = sum1 * sum2;
     }
 

        fixedScalingFactors[patIdx * PATTERN_BLOCK_SIZE + tx] =
            scalingFactors[blockIdx.x * PATTERN_BLOCK_SIZE * 4 + patIdx * PATTERN_BLOCK_SIZE + tx];

    if (patIdx == 0) {
        sMatrix1[tx] = matrix1[tx]; // All coalesced memory reads
        sMatrix2[tx] = matrix2[tx];
    }

                                                                REAL* matrices1,
                                                                REAL* matrices2,
                                                                int totalPatterns) {
    REAL sum1 = 1;
    REAL sum2;
    int i;

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