function I = flat_index(n,i,j)
    I = i + (j - 1)*n;
end