Hi Alex —
I didn't see your response earlier, but I reached the same conclusion about a half hour after you, and I can now explain what's going wrong. Happily (for me, as the implementor of the swap), it's not a bug in the swap, but rather a subtle issue in the user code—though I'm wondering whether perhaps we could've done more to protect you. It really is a subtle case.
Consider the following lines, running on a 10x10 array mapped to 2 locales, such that locale 0 owns {1..5, 1..10} and locale 1 owns {6..10, 1..10}:
for i in 1..rows {
T[i,cols+1] = 80.0*i/rows; // right-side boundary
Tnew[i,cols+1] = T[i,cols+1];
}
And keep in mind that with the fluff value of 1, each will actually own {0..6, 0..11} and {5..11, 0..11}, respectively.
Now, consider the iteration when i==5. Because this is a loop over a range (1..5) and ranges are not distributed, all iterations of the loop will be executed by the task that encountered the loop from the locale where it's executing. Since this loop doesn't appear within an on-clause, either explicitly or implicitly, it'll be executed on locale 0. I noticed that this would be a scalability problem as rows and cols get larger, but I didn't anticipate that it would also be a correctness problem (nor did I really think about these lines of code much at all, admittedly, since they seemed innocuous).
The problem is that the semantics of StencilDist are that if a value can be read locally, it will be. So what happens is that the read of T[i, cols+1] for i==5 is a read of T[6, 11], which falls within locale 0's fluff. Meanwhile, the previous write to T[i, cols+1] appropriately wrote to locale 1's copy of T[6, 11] since writes target the original copy of the array value; and, of course, we haven't called updateFluff() between the assignments (nor should we I think… feels like overkill). So the problem is that we're reading a default-initialized fluff value, getting the value 0.0 back, and storing that into Tnew[6, 11] on locale 1.
For the first iteration of the stencil, this isn't a problem since we're not reading Tnew anyway. For the versions that don't swap, it's also not a problem because we'll never read Tnew[6, 11]. However, for the version that swaps using the pointer-swap implementation, Tnew[6,11] becomes T[6,11] on the second iteration, and by the time we hit iteration 3, it starts resulting in errors propagating into the delta value. Specifically, if I print out the array on every iteration, I find that the rightmost computed column of values ceases to be monotonically increasing due to the 0.0 bringing its average down.
So, one solution would be to repeat the boundary condition expression in both assignments to avoid the read of the fluff, and this is why your constant boundary condition case also worked:
for i in 1..rows {
T[i,cols+1] = 80.0*i/rows; // right-side boundary
Tnew[i,cols+1] = 80.0*i/rows;
}
But a better approach for the sake of scalability would be to use a forall loop over a distributed domain or array to assign the boundaries since this will not only use multiple tasks, but also causes each task to compute that locale's local sub-array of values. Here's my preferred way to express that:
forall (i,j) in Dist.exterior(0,1) {
T[i,j] = 80.0*i/rows; // right-side boundary
Tnew[i,j] = T[i,j];
}
forall (i,j) in Dist.exterior(1,0) {
T[i,j] = 80.0*j/cols; // bottom-side boundary
Tnew[i,j] = T[i,j];
}
Trying to think of things we could've done to help alert you to this issue, I've been coming up short so far. I'd be curious if you, or others reading, have ideas.
I will be curious what kind of performance you get from this version once correct. And that reminds me that I also wanted to mention that my expectation (without a high degree of confidence) is that the previous approach to the reduction in which you compute it alongside the stencil seems like it ought to be faster at large problem sizes since it only requires traversing the pair of arrays once rather than twice, resulting in better locality.
Hope this all makes sense, and please ask questions if it doesn't.
-Brad