 #pragma unroll
         for (k = 0; k <= SHIFT; k++) {
           prods[j][SHIFT+k] += aa[j] * bb[j+k];
+        }
+#pragma unroll
+        for (k = 1; k <= SHIFT; k++) {
           prods[j][SHIFT-k] += aa[j+k] * bb[j];
         }
       }

      for (j = 0; j < WLEN; j++) {
#pragma unroll
        for (k = 0; k < window; k++) {
          CC[j + WLEN * threadIdx.y][k] = prods[j][k];
        }
      }
    }
    __syncthreads();
    for (j = 0; j < WLEN * window; j += dxy) {
      ic = (j + tid) / window;
      ir = (j + tid) - ic * window;
      for (k = 1; k < blockDim.y; k++) {
        if (j + tid < WLEN * window) {
          CC[ic][ir] += CC[ic + k * WLEN][ir];
        }
        __syncthreads();
      } 
      if (j + tid < WLEN * window && ic + jcol < ncols) {
        //        int x = ir + (ic + jcol) * window;
        //        if (x > window * ncols) printf("ir %d, ic %d, jcol %d\n", ir, ic, jcol);
        C[ir + (ic + jcol) * window] = CC[ic][ir];
      }
      }
  }
}


  float bb[WLEN + SHIFT];
  float prods[WLEN][window];
  int word[WLEN + SHIFT];
  __shared__ float CC[WLEN*BYDIM*window];
  int i, j, k, icol, jcol;
  int tid = threadIdx.x + blockDim.x * threadIdx.y;
  int dxy = blockDim.x * blockDim.y;

