Polytope Model

This is really cool and good to know about. I hadn't heard of this technique before, but have caught some interest in SPIRV lately while learning about WebGPU. Well I just noticed that Polly is used in the LLVM SPIRV backend... Makes sense! https://github.com/KhronosGroup/LLVM-SPIRV-Backend/blob/mast...

I'm very excited about the potential for writing GPU code from whatever high level language is convenient. (Which will not be WSL, by the way).

This is hard to grasp without a diagram. The detailed example is probably easier to understand.

The result computed in the current iteration is a[i - j][j]. The result computed in the iteration before that is a[i - (j - 1)][j - 1].

The example is super confusing because they define the alternative coordinate system (i', j') in the previous paragraph, but then proceed to use i and j (unprimed) in the code example, presumably because C does not allow for primes in variable names. p and q would have been better.

But, note that if you execute the inner loop for two consecutive values, j' and j' + 1, then you get the following computations:

a[i' - j'][j'] <- a[i' - j' - 1][j' ] + a[i' - j' ][j' - 1]

a[i' - j' - 1][j' + 1] <- a[i' - j' - 2][j' + 1] + a[i' - j' - 1][j']

You can see from the above that the left hand side of the first line does not actually appear on the right hand side of the second. This is because both indexes into a vary as j' changes.

It's a cool optimization, but I think the article fails to clearly explain the general method of how a compiler would arrive there.

Edit: I completely screwed up the indices in my initial post. They should be fixed now.

  In [1]: N = 5
  
  In [2]: def skewed():
     ...:     for p in range(2, 2*N - 1):
     ...:         tmin = max(1, p - (N - 1))
     ...:         tmax = min((p + 1) // 2, N - 1)
     ...:         for t in range(tmin, tmax + 1):
     ...:             print("p: {} t: {} | a[{}][{}] = a[{}][{}] + a[{}][{}]".format(
     ...:                 p,
     ...:                 t,
     ...:                 p - t,
     ...:                 t,
     ...:                 p - t - 1,
     ...:                 t,
     ...:                 p - t,
     ...:                 t - 1,
     ...:             ))
  
  In [3]: def unskewed():
     ...:     for i in range(1, N):
     ...:         for j in range(1, min(i+2, N)):
     ...:              print("i: {} j: {} | a[{}][{}] = a[{}][{}] + a[{}][{}]".format(
     ...:                  i,
     ...:                  j,
     ...:                  i,
     ...:                  j,
     ...:                  i - 1,
     ...:                  j,
     ...:                  i,
     ...:                  j - 1,
     ...:              ))
  
  In [4]: unskewed()
  i: 1 j: 1 | a[1][1] = a[0][1] + a[1][0]
  i: 1 j: 2 | a[1][2] = a[0][2] + a[1][1]
  i: 2 j: 1 | a[2][1] = a[1][1] + a[2][0]
  i: 2 j: 2 | a[2][2] = a[1][2] + a[2][1]
  i: 2 j: 3 | a[2][3] = a[1][3] + a[2][2]
  i: 3 j: 1 | a[3][1] = a[2][1] + a[3][0]
  i: 3 j: 2 | a[3][2] = a[2][2] + a[3][1]
  i: 3 j: 3 | a[3][3] = a[2][3] + a[3][2]
  i: 3 j: 4 | a[3][4] = a[2][4] + a[3][3]
  i: 4 j: 1 | a[4][1] = a[3][1] + a[4][0]
  i: 4 j: 2 | a[4][2] = a[3][2] + a[4][1]
  i: 4 j: 3 | a[4][3] = a[3][3] + a[4][2]
  i: 4 j: 4 | a[4][4] = a[3][4] + a[4][3]
  
  In [5]: skewed()
  p: 2 t: 1 | a[1][1] = a[0][1] + a[1][0]
  p: 3 t: 1 | a[2][1] = a[1][1] + a[2][0]
  p: 3 t: 2 | a[1][2] = a[0][2] + a[1][1]
  p: 4 t: 1 | a[3][1] = a[2][1] + a[3][0]
  p: 4 t: 2 | a[2][2] = a[1][2] + a[2][1]
  p: 5 t: 1 | a[4][1] = a[3][1] + a[4][0]
  p: 5 t: 2 | a[3][2] = a[2][2] + a[3][1]
  p: 5 t: 3 | a[2][3] = a[1][3] + a[2][2]
  p: 6 t: 2 | a[4][2] = a[3][2] + a[4][1]
  p: 6 t: 3 | a[3][3] = a[2][3] + a[3][2]
  p: 7 t: 3 | a[4][3] = a[3][3] + a[4][2]
  p: 7 t: 4 | a[3][4] = a[2][4] + a[3][3]
  p: 8 t: 4 | a[4][4] = a[3][4] + a[4][3]

Actually writing it out like that really helped me grok what was going on, personally. If you look at the unskewed function then you can see that for every i there is a value for j that relies on the output from a previous j, e.g. the i: 1 j: 2 iteration uses a[1][1] which is the output from the i: 1 j: 1 iteration. But you can't find a similar example in the skewed function.