I need to solve a linear system $Ax=b$ many times (billions), each time with small changes in A and b from the previous time. The solution can be approximate, but it needs to be very fast: if $x$ is of size $N$ (typically $500$ to $1000$), there is only time to do around $100 N^2$ floating point operations.
The matrix A has a special structure: it is symmetrical, and the rows and columns correspond to points embedded in 3D; the off-diagonal values are $c_{ij}/r_{ij}$ where $r_{ij}$ is the distance and $c_{ij}$ is some coefficient (the coefficients are not the same but they are all within a relatively narrow range, let's say +-1). The density of points in 3D (number of points per unit volume) is also roughly constant. The changes in A correspond to small motions of the points in 3D space: only $r_{ij}$ changes, $c_{ij}$ stays constant.
It seems that the points in 3D are really important for an efficient solution to this. If a single point moves, the x values of the solution for it and its nearby neighbors in 3D (which have large matrix elements) would change significantly, but for distant neighbors the change would be very small, both because the matrix elements on the order of $1/r$ are small and because the relative change in those elements is small. It should be possible to come up with an efficient algorithm which takes iterative steps towards an (approximate) solution but essentially propagates changes based on distance; so in calculating $Ax$, for each $x_i$ the $\sum A_{ij} x_{j}$ for distant j is updated a lot less often (can be stored and reused for many steps) while for nearest neighbor j it is updated every step.
Is there already a known algorithm for something like this?
P.S. This seems closely related to N-body problems, esp tree methods - by analogy a saved $\sum A_{ij} x_{j}$ is equivalent to the single particle which represents a distant branch of the tree in Barnes-Hut; and also to hierarchical time step methods - by analogy updating the saved value less often corresponds to using a longer time step for distant interactions. The main difference seems to be that I need to find an approximate solution to a linear system (ideally with some guarantees on accuracy and/or convergence to the correct solution if the particles stop moving), not take an integrator time step.
P.S. Probably worth a mention, while the number of arithmetic operations to find an approximate solution is limited, the amount of additional data stored together with the previous solution can be very large; much larger than N^2 or even N^3