   float *weight_data = THCudaTensor_data(kernel);
   float *output_data = THCudaTensor_data(output);
 
-  // auto compute nb of blocks and threads
+  // auto compute nb of blocks and threads:
+  // this might look a bit arbitrary, but it's not. We want to be in a sweet spot,
+  // where most of the time we'll have 16x16 threads per output map, that is 256 threads
+  // dealing with each map, and one block per map.
+  // example: if we have 16 output maps, each being 256x256, we'll have 16 blocks
+  // with 256 threads, each thread processing a 16x16 suboutput.
+  // call me if it's not clear ;-)
+  int patch_w = (int)(pow(2, ceil(log2((float)nOutputCols))) / 16);
+  if (patch_w < 2) patch_w = 2;
+  else if (patch_w > 32) patch_w = 32;
+  int patch_h = patch_w;
   dim3 blocks(nOutputPlane);
-  dim3 threads(1, nOutputRows / 4);
-  if ((nOutputRows % 4) != 0) threads.y++;
+  dim3 threads(nOutputCols / patch_w, nOutputRows / patch_h);
+  if ((nOutputRows % patch_h) != 0) threads.y++;
+  if ((nOutputCols % patch_w) != 0) threads.x++;
 
   // convolution: input with kernel, 4 modes (full, valid, xcorr2 or conv2)
   if (type[0] == 'f') {

  } else { // 'v'

    if (type[1] == 'x') {
      validXCorr2ptr <<<blocks, threads>>> (input_data, weight_data, output_data,
                                            nInputPlane, nInputRows, nInputCols,
                                            nOutputPlane*nInputPlane, nKernelRows, nKernelCols,
                                            srow, scol);
    } else {
      THError("valid conv2 not implemented yet");
    }

  }

#include "THCGeneral.h"

/*
 * Description:
 *   This code provides convolutions and xcorrelations that are API compatible with
 *   the ones in THLabConv.
 *
 * History:
 *   July 21, 2011, 11:21PM  -  Clement Farabet  -  Creation, based conv2d routine
 */

#define CUDA_SHARED_MEM_SIZE (4*1024-32) // this is given by nVidia: max shared mem per block

/*
 * Description:
 *   base conv2D routine: 3D input, 3D output, 4D kernel
 *
 *   - all chunks of data should be contiguous
 *   - the swapkernel flag can be used to generate a conv2 instead of xcorr2
 *   - the templated kernel size is useful to generate code that's 2x faster
 *     but can be set to 0 to allow arbitrary kernel sizes
 *
 * to do:
 *   + the most crucial thing is to make the threads make contiguous
 *     memory accesses wrt each other, to enable memory reads to coalesce.
 *     this should give an extra 4 to 8x speedup
 */
template <bool swapkernel, int T_kernel_h, int T_kernel_w>
  __global__ void conv2generic(float *input, float *kernel, float *output,
                               int input_n, int input_h, int input_w,
                               int kernel_n, int kernel_h, int kernel_w,
                               int stride_h, int stride_w,
                               int patch_h, int patch_w)
{
  // output dimensions
  int output_h = (input_h - kernel_h) / stride_h + 1;
  int output_w = (input_w - kernel_w) / stride_w + 1;

  // xcorr or conv
  int koffset = swapkernel ? kernel_w*kernel_h-1 : 0;

  // generate offsets according to block/thread ids
  int xx_start = threadIdx.x * patch_w;
  int xx_end = xx_start + patch_w; if (xx_end > output_w) xx_end = output_w;
  int yy_start = threadIdx.y * patch_h;
  int yy_end = yy_start + patch_h; if (yy_end > output_h) yy_end = output_h;
  int oo_start = blockIdx.x;
  int oo_end = oo_start+1;
  int ii_start = 0;
  int ii_end = input_n;

  // iterators
  int oo, ii, xx, yy, kx, ky, kk;

  // do the kernels fit in shared mem ?
  if (input_n*kernel_w*kernel_h <= CUDA_SHARED_MEM_SIZE) {

    // put the kernel in shared memory
    __shared__ float shared_kernel[CUDA_SHARED_MEM_SIZE];

    // first thread of each block does the copy
    if ((threadIdx.x == 0) && (threadIdx.y == 0)) {
      for (kk = 0; kk < kernel_w*kernel_h*input_n; kk++) {
        shared_kernel[kk] = kernel[input_n*kernel_w*kernel_h*blockIdx.x + kk];
      }
    }

    // sync threads
    __syncthreads();

    // templated kernel size
    if ((T_kernel_w > 0) && (T_kernel_h > 0)) {
      // unrolled convolution loop
      for(oo = oo_start; oo < oo_end; oo++) {
        for(ii = ii_start; ii < ii_end; ii++) {
          for(yy = yy_start; yy < yy_end; yy++) {
            for(xx = xx_start; xx < xx_end; xx++) {
              // Dot product in two dimensions... (between input image and the mask)
              float *input_p = input + ii*input_h*input_w + yy*stride_h*input_w + xx*stride_w;
              float *output_p = output + oo*output_h*output_w + yy*output_w + xx;
              float *kernel_p = shared_kernel + ii * kernel_w * kernel_h + koffset;
              float sum = 0;
              if (swapkernel) {
#pragma unroll
                for(ky = 0; ky < T_kernel_h; ky++) {
#pragma unroll
                  for(kx = 0; kx < T_kernel_w; kx++) {
                    sum += input_p[kx]*(*kernel_p--);
                  }
                  input_p += input_w;
                }
              } else {
#pragma unroll
                for(ky = 0; ky < T_kernel_h; ky++) {
#pragma unroll
                  for(kx = 0; kx < T_kernel_w; kx++) {
                    sum += input_p[kx]*(*kernel_p++);
                  }
                  input_p += input_w;
                }
              }
              *output_p += sum;
            }
          }
        }
      }
    } else {
      // default convolution loop
      for(oo = oo_start; oo < oo_end; oo++) {
        for(ii = ii_start; ii < ii_end; ii++) {
          for(yy = yy_start; yy < yy_end; yy++) {
            for(xx = xx_start; xx < xx_end; xx++) {
              // Dot product in two dimensions... (between input image and the mask)
              float *input_p = input + ii*input_h*input_w + yy*stride_h*input_w + xx*stride_w;
              float *output_p = output + oo*output_h*output_w + yy*output_w + xx;
              float *kernel_p = shared_kernel + ii * kernel_w * kernel_h + koffset;
              float sum = 0;
              for(ky = 0; ky < kernel_h; ky++) {
                for(kx = 0; kx < kernel_w; kx++) {
                  if (swapkernel) sum += input_p[kx]*(*kernel_p--);
                  else sum += input_p[kx]*(*kernel_p++);
                }
                input_p += input_w;
              }
              *output_p += sum;
            }
          }
        }
      }
    }

  } else { // not enough shared mem for kernels, simply stream them

    // convolution loop
    for(oo = oo_start; oo < oo_end; oo++) {
      for(ii = ii_start; ii < ii_end; ii++) {
        for(yy = yy_start; yy < yy_end; yy++) {
          for(xx = xx_start; xx < xx_end; xx++) {
            // Dot product in two dimensions... (between input image and the mask)
            float *input_p = input + ii*input_h*input_w + yy*stride_h*input_w + xx*stride_w;
            float *output_p = output + oo*output_h*output_w + yy*output_w + xx;
            float *kernel_p = kernel + (oo * input_n + ii) * kernel_w * kernel_h + koffset;
            float sum = 0;
            for(ky = 0; ky < kernel_h; ky++) {
              for(kx = 0; kx < kernel_w; kx++) {
                if (swapkernel) sum += input_p[kx]*(*kernel_p--);
                else sum += input_p[kx]*(*kernel_p++);
              }
              input_p += input_w;
            }
            *output_p += sum;
          }
        }
      }
    }
  }
}

/*
 * API-compatible with THRealTensor_conv2Dmv
 * 3D input, 4D kernel, 3D output
 * matrix vector product like: y <- Ax + beta*y
 */
TH_API void THCudaTensor_conv2Dmv(THCudaTensor *output, float beta, THCudaTensor *input,
                                  THCudaTensor *kernel, long srow, long scol, const char *type)
{
