Fix the indexing #101

The old index x[i1 + d1 - s_in] only works when the matrix has 2 * p_in + 1 diagonals. Diagonal d1 of row i1 is the column i1 - q + d1, where q is the pad of the matrix. That column sits at x[i1 - q + d1 - s_in + p_in], because x is padded with p_in. When q = p_in this reduces to the old i1 + d1 - s_in. When the matrix has fewer diagonals (q < p_in, e.g. a block between spaces of different degree), the old index is shifted by p_in - q, which gives wrong results. That's why those matrices used to raise NotImplementedError.

So off1 = p_in - q - s_in is that shift. It is computed once outside the loop, and for q = p_in it is exactly -s_in, the old offset.

I agree off1 hides what it means. I can write it out instead, e.g.

q1 = (mat.shape[1] - 1) // 2  # pads of the matrix
...
val += mat[p_out + i1_loc, d1] * x[i1 - q1 + d1 - s_in + p_in]

so the column i1 - q1 + d1 can be read directly from the index. Would that be clearer?

Originally posted by @max-models in #97 (comment)