   const unsigned nThreadPerBlock = 256;
   unsigned nBlockPerColumn = (size[0] + nThreadPerBlock - 1) / nThreadPerBlock;
   dim3 threads(nThreadPerBlock);
-  dim3 grid(nBlockPerColumn, size[1], size[2]);
+  unsigned maxGridDim = 1024; // anything < 64k is fine. The choice has no impact on performance.
+  dim3 grid(min(maxGridDim, nBlockPerColumn), min(maxGridDim, size[1]), min(maxGridDim, size[2]));
 
   THCudaTensor_kernel_reduceOuterDim<<<grid, threads>>>(THCudaTensor_data(tgt),
           THCudaTensor_data(src), src_stride, tgt_stride, size, op, init);

        dim4 src_stride, dim4 tgt_stride, dim4 size, Op op, T init)
{
  const size_t reduce = 3;
  unsigned col = blockIdx.x * blockDim.x + threadIdx.x;
  float *src = src_ + blockIdx.z * src_stride[2] + blockIdx.y * src_stride[1] + col;

  if(col < size[0]) {
    T acc = init;
    for(unsigned i=0; i < size[reduce]; i++) {
      acc = op((T) acc, (T) *src);
      src += src_stride[reduce];
    }
    tgt[blockIdx.z * tgt_stride[2] + blockIdx.y * tgt_stride[1] + col] = float(acc);
  }
}


template<class Op, typename T>


  dim3 threads(32, 16);
  unsigned nBlockPerRow = (size[1] + threads.y - 1) / threads.y;
  unsigned maxGridDim = 1024; // anything < 64k is fine. The choice has no impact on performance.
  dim3 grid(min(maxGridDim, size[2]), min(maxGridDim, nBlockPerRow), min(maxGridDim, size[3]));


  THCudaTensor_kernel_reduceInnermostDim<<<grid, threads>>>(THCudaTensor_data(tgt),
          THCudaTensor_data(src), src_stride, tgt_stride, size, op, init);
