 
 // -----------------------------------------------------------------------------
 
+#define USE_NVIDIA_OPTIMIZATION
+#ifdef USE_NVIDIA_OPTIMIZATION
+
+template< int NUM_ELEMENTS, int NUM_THREADS_PER_BLOCK >
+__global__ void computeStencilsNv(float const *__restrict cvs,
+                                  float * vbuffer,
+                                  unsigned char const *__restrict sizes,
+                                  int const *__restrict offsets,
+                                  int const *__restrict indices,
+                                  float const *__restrict weights,
+                                  int start,
+                                  int end)
+{
+  // Shared memory to stage indices/weights.
+  __shared__ int   smem_indices_buffer[NUM_THREADS_PER_BLOCK];
+  __shared__ float smem_weights_buffer[NUM_THREADS_PER_BLOCK];
+
+  // The size of a single warp.
+  const int WARP_SIZE = 32;
+  // The number of warps per block.
+  const int NUM_WARPS_PER_BLOCK = NUM_THREADS_PER_BLOCK / WARP_SIZE;
+  // The number of outputs computed by a single warp.
+  const int NUM_OUTPUTS_PER_WARP = WARP_SIZE / NUM_ELEMENTS;
+  // The number of outputs computed by a block of threads.
+  const int NUM_OUTPUTS_PER_BLOCK = NUM_OUTPUTS_PER_WARP*NUM_WARPS_PER_BLOCK;
+  // The number of active threads in a warp.
+  const int NUM_ACTIVE_THREADS_PER_WARP = NUM_OUTPUTS_PER_WARP * NUM_ELEMENTS;
+
+  // The number of the warp inside the block.
+  const int warpId = threadIdx.x / WARP_SIZE;
+  const int laneId = threadIdx.x % WARP_SIZE;
+
+  // We use NUM_ELEMENTS threads per output. Find which output/element a thread works on.
+  int outputIdx = warpId*NUM_OUTPUTS_PER_WARP + laneId/NUM_ELEMENTS, elementIdx = laneId%NUM_ELEMENTS;
+
+  // Each output corresponds to a section of shared memory.
+  volatile int   *smem_indices = &smem_indices_buffer[warpId*WARP_SIZE + (laneId/NUM_ELEMENTS)*NUM_ELEMENTS];
+  volatile float *smem_weights = &smem_weights_buffer[warpId*WARP_SIZE + (laneId/NUM_ELEMENTS)*NUM_ELEMENTS];
+
+  // Disable threads that have nothing to do inside the warp.
+  int i = end;
+  if( laneId < NUM_ACTIVE_THREADS_PER_WARP )
+    i = start + blockIdx.x*NUM_OUTPUTS_PER_BLOCK + outputIdx;
+
+  // Iterate over the vertices.
+  for( ; i < end ; i += gridDim.x*NUM_OUTPUTS_PER_BLOCK )
+  {
+    // Each thread computes an element of the final vertex.
+    float x[NUM_ELEMENTS == 4 ? 4 : 1] = {0.f};
+
+    // Load the offset and the size for each vertex. We have NUM_THREADS_PER_VERTEX threads loading the same value.
+    const int offset_i = offsets[i], size_i = sizes[i];
+
+    // Iterate over the stencil.
+    for( int j = offset_i, j_end = offset_i+size_i ; j < j_end ; )
+    {
+      int j_it = j + elementIdx;
+
+      // Load some indices and some weights. The transaction is coalesced.
+      smem_indices[elementIdx] = j_it < j_end ? indices[j_it] : 0;
+      smem_weights[elementIdx] = j_it < j_end ? weights[j_it] : 0.f;
+
+      // Thread now collaborates to load the vertices.
+      if( NUM_ELEMENTS == 4 )
+      {
+        if( j_it < j_end )
+        {
+          float w = smem_weights[elementIdx];
+          float4 tmp = reinterpret_cast<const float4 *__restrict>(cvs)[smem_indices[elementIdx]];
+          x[0] += w*tmp.x;
+          x[1] += w*tmp.y;
+          x[2] += w*tmp.z;
+          x[3] += w*tmp.w;
+        }
+        j += 4;
+      }
+      else
+      {
+        #pragma unroll
+        for( int k = 0 ; k < NUM_ELEMENTS ; ++k, ++j )
+          if( j < j_end )
+            x[0] += smem_weights[k] * cvs[smem_indices[k]*NUM_ELEMENTS + elementIdx];
+      }
+    }
+
+    // Store the vertex.
+    if( NUM_ELEMENTS == 4 )
+      reinterpret_cast<float4*>(vbuffer)[i] = make_float4(x[0], x[1], x[2], x[3]);
+    else
+      vbuffer[NUM_ELEMENTS*i + elementIdx] = x[0];
+  }
+}
+
+#endif USE_NVIDIA_OPTIMIZATION
+
+// -----------------------------------------------------------------------------
+
 #include "../version.h"
 
 #define OPT_KERNEL(NUM_ELEMENTS, KERNEL, X, Y, ARG) \

    for(int i = 0; i < count; ++i) dst[i] += src[i] * weight;
}

// -----------------------------------------------------------------------------

template <int NUM_ELEMENTS> __global__ void
computeStencils(float const * cvs, float * vbuffer,

        return;
    }

#ifdef USE_NVIDIA_OPTIMIZATION
    OPT_KERNEL_NVIDIA(3, computeStencilsNv, 2048, 256, (cvs, dst, sizes, offsets, indices, weights, start, end));
    OPT_KERNEL_NVIDIA(4, computeStencilsNv, 2048, 256, (cvs, dst, sizes, offsets, indices, weights, start, end));
#else
    OPT_KERNEL(3, computeStencils, 512, 32, (cvs, dst, sizes, offsets, indices, weights, start, end));
    OPT_KERNEL(4, computeStencils, 512, 32, (cvs, dst, sizes, offsets, indices, weights, start, end));
#endif

    computeStencils <<<512, 32>>>(cvs, dst, length, stride,
        sizes, offsets, indices, weights, start, end);
