   int i, ioff, ioffmv, ip, j, k, c, tmp, cmin, imin;
   int zero = 0;
   int sid = threadIdx.x + blockDim.x * threadIdx.y;
+
+  // Load data into shared memory
   for (i = 0; i < TLEN/1024; i++) {
     sa[sid * i*1024] = a[sid + i*1024 + TLEN*blockIdx.x];
   }

      imin = __shfl_down(imin, k);
    }
  }
  if (threadIdx.x == 0) {
    sop[threadIdx.y] = imin;
    sow[threadIdx.y] = cmin;
  }
  __syncthreads();
  if (threadIdx.y == 0) {
    op[threadIdx.x + 32*blockIdx.x] = sop[threadIdx.x];
    ow[threadIdx.x + 32*blockIdx.x] = sow[threadIdx.x];

    bb[j] = sb[threadIdx.y][j + ioffmv];
    ww[j] = sw[threadIdx.y][j + ioffmv];
  }

  // Step through offsets in A string
  for (j = 0; j < TLEN*NVEC/8; j++) {
    tmp = VECLEN + ioff + j / 4;
    if (tmp - ioffmv < TLEN - VECLEN * NVEC) {
      if (j % 4 == 0) {
        aa[VECLEN] = sa[tmp];
      }
      c = 0;
      // Inner loop over the length of the vector in registers
#pragma unroll
      for (k = 0; k < VECLEN; k++) {    
        hammingcell(aa[k], aa[k+1], bb[k], ww[k], c, tmp, zero);
      }
      hammingend(aa[k]);
      // Need to sum over NVEC to get complete score for a string
#pragma unroll
      for (k = 1; k < NVEC; k *= 2) {    
        c = c + __shfl_down(c, k);
      }
      // Now compare with the accumulated min
      if (c < cmin) {
        cmin = c;
        imin = 4 * ioff + j;
      }
    }
  }
  // Compute the min across groups of NVEC threads in this warp
  for (k = NVEC; k < 32; k *= 2) {    
    tmp = __shfl_down(cmin, k);
    if (tmp < cmin) {
