   );
 
   THCudaCheck(cudaFree(stride_));
+  THCudaTensor_free(state, indices);
+}
+
+__global__ void THCudaTensor_kernel_indexSelect_contiguous(
+  float *tensor, float *src, long stride, float *index, long idxSize)
+{
+  // In the typical case, each block of 128 threads handles a 4x128
+  // section of the output with each warp handling a single 1x128 row.
+  // The outer loops handle inputs larger than 4*65535 or strides larger
+  // than 128*65535.
+  const int VT = 4;
+  const int WARP_SIZE = 32;
+  const int MAX_DIM_SIZE = 65535;
+
+  for (int idx = blockIdx.x * blockDim.y + threadIdx.y; idx < idxSize; idx += blockDim.y * MAX_DIM_SIZE) {
+    for (int startIdx = threadIdx.x + blockIdx.y * VT*WARP_SIZE; startIdx < stride; startIdx += VT*WARP_SIZE*MAX_DIM_SIZE) {
+      const int srcIdx = ((int) index[idx] - 1) * stride;
+      const int targetIdx = idx * stride;
+
+      #pragma unroll
+      for (int i = 0; i < VT; i++) {
+        const int featureIdx = startIdx + i * WARP_SIZE;
+        if (featureIdx < stride) {
+          tensor[targetIdx + featureIdx] = src[srcIdx + featureIdx];
+        }
+      }
+    }
+  }
 }
 
 __global__ void THCudaTensor_kernel_indexSelect(

  THArgCheck(nIndex == src->size[dim], 4, "length of src.size[dim] is not equal to length of indices");

  src = THCudaTensor_newContiguous(state, src);

  nRes = THCudaTensor_nElement(state, res_);
  dim3 nthreads(16, 16);

  THCudaTensor_resize(state, res_, newSize, NULL);
  THLongStorage_free(newSize);

  res = THCudaTensor_newContiguous(state, res_);
  indices = THCudaTensor_newContiguous(state, indices);

  if (THCudaTensor_isContiguous(state, src) && dim == 0)
  {
    long stride = src->stride[0];

    int blockX = std::min(DIVUP(nIndex, 4), 65535L);
    int blockY = std::min(DIVUP(stride, 128), 65535L);

    dim3 nthreads(32, 4);
    dim3 nblocks(blockX, blockY);

    THCudaTensor_kernel_indexSelect_contiguous<<<nblocks, nthreads, 0, stream>>>(
      THCudaTensor_data(state, res),
      THCudaTensor_data(state, src),
      stride,
      THCudaTensor_data(state, indices),
      nIndex);

    THCudaTensor_free(state, indices);
    THCudaTensor_freeCopyTo(state, res, res_);
    return;
  }

  nRes = THCudaTensor_nElement(state, res);
  dim3 nthreads(16, 16);
  dim3 nblocks(ceil((float)nRes / nIndex / (16*16)));

  THCudaCheck(cudaMalloc((void**)&stride_, src->nDimension * sizeof(long)));
  THCudaCheck(cudaMemcpy(stride_, src->stride, src->nDimension * sizeof(long), cudaMemcpyHostToDevice));

  THCudaTensor_kernel_indexSelect<<<nblocks, nthreads, 0, stream>>>(
    THCudaTensor_data(state, res), THCudaTensor_data(state, src),
    stride_, THCudaTensor_data(state, indices),
    src->nDimension, dim, nIndex, nRes, src->size[dim]
  );

  THCudaCheck(cudaFree(stride_));
  THCudaTensor_free(state, indices);
  THCudaTensor_freeCopyTo(state, res, res_);
}
