     /* BPROP                                                                 */
     /*************************************************************************/
 
-    /* How many output layers? */
+    /* Make this_layer_outputs point to the outputs of the layer before the
+       output */
     this_layer_outputs -= no;
 
     layer_weights -= ni * w_stride;
 
-    /* First error calculation pass */
-    float last_output[N];
-    FOR_ALL_X(last_output[x] = scratch[x * scratch_stride + tid]);
+    /* We keep the outputs of the next layer here to avoid needing to reload
+       them. */
+    float prev_outputs[N];
+    FOR_ALL_X(prev_outputs[x] = scratch[x * scratch_stride + tid]);
 
     __syncthreads();
 
+    /* Calculate the error terms for the output units.  Scratch contains the
+       outputs of output layer, so we can access this instead of
+       this_layer_outputs. */
     for (int x = 0;  x < N && x < valid_examples;  ++x) {
         bool correct = (labels[x] == tid);
         float wanted = (correct ? fire : inhibit);
         scratch[x * scratch_stride + tid]
-            = (tid < no ? wanted - last_output[x] : 0.0);
+            = (tid < no ? wanted - prev_outputs[x] : 0.0);
     }
     
-    /* Let everything catch up */
+    /* Let everything catch up.  Scratch now contains the errors for the
+       output layer. */
     __syncthreads();
 
     /* Backpropegate. */

        /* Now for the updates.  In order to avoid trying to write the same
           memory over and over, we stagger the starting points so that
           each example will start at a different place, thus minimising
           conflicting writes when we have multiple multiprocessors working
           on the same thing. */

        int thread_stride = ni / num_threads_in_block;
        if (thread_stride == 0) thread_stride = 1;

            double total_update = 0.0;

            for (int x = 0;  x < N && x < valid_examples;  ++x) {
                // All threads read the same value
                float prev
                    = (l == 0
                       ? input[x * input_stride + i]
                       : last_layer_outputs[x * total_neurons + i]); 
                float update = prev * k[x] * d[x];
                total_update += update;

                // If we've just read the value we need for the next
                // block, then store it.
                if (i == tid) prev_outputs[x] = prev;
            }
            
            atomic_add(layer_updates[i * w_stride + tid], total_update);
