Consider the following reweighted Poisson equation: given $\operatorname{Q} $ and $g$, $$ (\nabla \cdot \operatorname{Q} \nabla ) f = g, $$ where $\operatorname{Q} $ is a diagonal matrix with positive entries. Note that $\nabla \cdot \nabla$ is the normal Laplace operator.
I am working on a discrete reweighted Poisson equation in high dimensions like 10000*10000. Since $\nabla \cdot \operatorname{Q} \nabla $ is a sparse matrix, I can use sparse linear solvers(the fastest package I found is scipy.spsolve) to solve it. But it's still not efficient enough. Regular FFT cannot be applied directly due to the existence of the reweighting matrix $\operatorname{Q}$. I wonder if anyone knows a better way to solve it, possibly other adaptions of FFT or specific sparse linear solvers?