@smth I believe the minimal example provided in my original post is the script you are asking about.  You just have to remove all the .cuda() if you want a CPU version.  That said, I’ve been playing around with that code this morning, and I’m finding all sorts of weird behavior.
First, the most serious problem is that nn.MaxPool3d outputs junk indices.  Here is an example:
pool3d = nn.MaxPool3d(kernel_size=2,stride=2,return_indices=True)
img3d = Variable(torch.rand(1,1,4,4,4))
out, indices = pool3d(img3d)
print(indices)

The elements of indices should be in the range [0,63], but here is what is printed:
Variable containing:
(0 ,0 ,0 ,.,.) = 
  4.6117e+18 -9.2234e+18
  3.0065e+10  4.2950e+09

(0 ,0 ,1 ,.,.) = 
  6.5536e+04  8.5213e+14
 -3.4588e+18  4.6117e+18
[torch.LongTensor of size 1x1x2x2x2]

Second, and this is the really “weird behavior” stuff, my minimal example stops producing an error when I change these lines of code:
input_batch = Variable(input_batch.cuda())
target_batch = Variable(target_batch.cuda())

to this:
input_batch = Variable(input_batch)
target_batch = Variable(target_batch)
input_batch = input_batch.cuda()
target_batch = target_batch.cuda()

I’m not sure what to make of that.  When I run the code on the CPU, I get an error complaining about the pooling indices, which makes sense.  However, when running on the GPU, I get no error at all, as long as I move the inputs and targets to the GPU in a very particular way (otherwise I get an illegal memory access error).  That’s dangerous behavior, because I wouldn’t have even noticed that the pooling indices are all messed up if I had written my code in a slightly different way.  (Another seemingly random thing is that removing the shortcut connection also eliminates the error, so this is all very weird.)
Finally, for future reference, regarding the correctness of my minimal example, I need to be more careful when slicing the pooling indices.  The pooling indices are supposed to refer to coordinates in the input Tensor given in flattened spatial coordinates.  So, what I really should have written is:
pre_pool_shape = out.size[2:]
out, indices = self.pool(out)
indices = indicies[:,:,pre_pool_shape[0]//2]   # only want the center slice
indices %= pre_pool_shape[1]*pre_pool_shape[2]  # flattened coords for 2d slice

Of course this won’t work until nn.MaxPool3d is fixed.