         const float * layer_weights = w[l];
         int w_stride = w_strides[l];
 
+        float * layer_updates = w_updates[l];
+        float * layer_bias_updates  = b_updates[l];
+        
         next_layer_outputs = this_layer_outputs - no;
         
         if (tid >= no) continue;
-
+        
         float d = delta(last_layer_outputs[tid], errors[tid], activation);
 
-        if (l > 1) {
-            /* Calculate new errors (for the next layer). */
+        /* Calculate the new error terms for the next layer */
+        // TODO: atomic... and then find a way to avoid data dependencies...
+        if (l > 1)
             for (unsigned i = 0;  i < ni;  ++i)
-                errors[i] = SIMD::vec_dotprod_dp(&delta[0],
-                                                 &layer.weights[i][0],
-                                                 no);
+                errors[i] += d * layer_weights[i * w_stride + tid];
+
+
+        /* Update the weights. */
+        float k = example_weight * learning_rate;
+        for (unsigned i = 0;  i < ni;  ++i) {
+            // No bank conflicts here as all threads are reading with the same
+            // i value
+            float k2 = last_layer_outputs[i] * k;
+            
+            layer_updates[i * w_stride + tid] += k2 * d;
         }
         
-
-        for (unsigned i = 0;  i < ni;  ++i) {
-            float inval = (l == 0 ? input[i] : last_layer_outputs[i]);
-            
-            // Coalesced access; maybe texture would be better
-            float weight = layer_weights[i * w_stride + tid];
-            
-            accum += weight * inval;
-        }
-
-        /* Update the weights. */
-        float k = w * learning_rate;
-        for (unsigned i = 0;  i < ni;  ++i) {
-            float k2 = layer_outputs[l - 1][i] * k;
-            SIMD::vec_add(&sub_weight_updates[l][i][0], k2, &delta[0],
-                          &sub_weight_updates[l][i][0], no);
-        }
-        
-        /* Update the bias terms.  The previous layer output (input) is
-           always 1. */
-        SIMD::vec_add(&sub_bias_updates[l][0], k, &delta[0],
-                      &sub_bias_updates[l][0], no);
-    }
-    
-#endif
-
+        /* Update the bias */
+        layer_bias_updates[tid] += k * d;
+    }
 }
+
+namespace ML {
+namespace CUDA {
+
+void train_examples(const float * feature_vectors,
+                    int num_feature_vectors,
+                    int feature_vector_width,
+                    const float * example_weights,
+                    const int * labels,
+                    int num_layers,
+                    const int * architecture,
+                    const float * const * weights,
+                    const float * const * biases,
+                    const int * w_strides,
+                    float * const * weight_updates,
+                    float * const * bias_updates,
+                    Activation activation,
+                    float fire,
+                    float inhibit,
+                    float learning_rate)
+{
+    DeviceData<float> d_feature_vectors
+        (feature_vectors,
+         num_feature_vectors * feature_vector_width);
+
+    if (feature_vector_width != architecture[0])
+        throw Exception("number of inputs doesn't match");
+
+    DeviceData<float> d_example_weights(example_weights,
+                                              num_feature_vectors);
+
+    DeviceData<int> d_labels(labels, num_feature_vectors);
+
+    DeviceData<int> d_architecture(architecture, num_layers + 1);
+
+    vector<DeviceData<float> > d_weights_storage(num_layers);
+    const float * d_weights[num_layers];
+
+    for (unsigned l = 0;  l < num_layers;  ++l) {
+        int no = architecture[l + 1];
+        int w_stride = w_strides[l];
+        d_weights_storage[l].init(weights[l], no * w_stride);
+        d_weights[l] = d_weights_storage[l];
+    }
+    
+    vector<DeviceData<float> > d_biases_storage(num_layers);
+    const float * d_biases[num_layers];
+
+    for (unsigned l = 0;  l < num_layers;  ++l) {
+        int no = architecture[l + 1];
+        d_biases_storage[l].init(biases[l], no);
+        d_biases[l] = d_biases_storage[l];
+    }
+    
+    DeviceData<int> d_w_strides(w_strides, num_layers);
+
+    vector<DeviceData<float> > d_weight_updates_storage(num_layers);
+    float * d_weight_updates[num_layers];
+
+    for (unsigned l = 0;  l < num_layers;  ++l) {
+        int no = architecture[l + 1];
+        int w_stride = w_strides[l];
+        d_weight_updates_storage[l].init(weights[l], no * w_stride);
+        d_weight_updates[l] = d_weight_updates_storage[l];
+    }
+
+    vector<DeviceData<float> > d_bias_updates_storage(num_layers);
+    float * d_bias_updates[num_layers];
+
+    for (unsigned l = 0;  l < num_layers;  ++l) {
+        int no = architecture[l + 1];
+        d_bias_updates_storage[l].init(biases[l], no);
+        d_bias_updates[l] = d_bias_updates_storage[l];
+    }
+
+    int max_width = 0;
+    int total_neurons = 0;
+    for (unsigned l = 0;  l <= num_layers;  ++l) {
+        max_width = max(max_width, architecture[l]);
+        total_neurons += architecture[l];
+    }
+
+    // We need our grid size to be exactly the maximum width of the output
+    dim3 threads(max_width);
+
+    // Our grid size is one per example
+    dim3 grid(num_feature_vectors);
+
+    size_t shared_mem_size = (max_width + total_neurons) * sizeof(float);
+
+    train_example_kernel <<< grid, threads, shared_mem_size >>>
+        (d_feature_vectors,
+         feature_vector_width,
+         d_labels,
+         d_example_weights,
+         num_layers,
+         d_weights,
+         d_biases,
+         d_architecture,
+         d_w_strides,
+         d_weight_updates,
+         d_bias_updates,
+         activation,
+         fire,
+         inhibit,
+         learning_rate);
+}
+
+} // namespace CUDA
+} // namespace ML

    __syncthreads();

    /* Backpropegate. */
    for (int l = nl - 1;  l >= 1;
         --l,
             __syncthreads(),
             last_layer_outputs = this_layer_outputs,

}

/* Given an output and an error, what's the delta? */
__device__ float delta(float output, float error, int activation)
{
    switch (activation) {
    case ML::ACT_TANH:
