                                                                       REAL* matrices2,
                                                                       REAL* scalingFactors,
                                                                       int totalPatterns) {
-    REAL sum1 = 0;
-    REAL sum2 = 0;
+    REAL sum1;
+    REAL sum2;
     int i;
 
     int tx = threadIdx.x;
-    int state = tx % 4;
+    int state = tx & 0x3; //  tx % 4;
     int pat = tx / 4;
     int patIdx = threadIdx.y;
     int matrix = blockIdx.y;

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

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

    __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];
//        }

//        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];
    }

