function y = idx(i,j,nj)
  y = (i-1)*nj + (j-1) + 1;
end