How to solve a problem with structure similar to a finite difference discretization of the 2D Poisson equation, but with non-symetric coefficients?

Recently, I've been asking about methods to solve a finite difference discretization of the 2D Poisson equation (see here and here) of the form:

$$U_{i-1,j} + U_{i+1,j} -4U_{i,j} + U_{i,j-1} + U_{i,j+1} = f_{i,j}$$

The left hand side of this system is often referred to as the 5-point discrete Laplacian. Suppose that left hand side of the linear system is not quite the traditional five-point discrete Laplacian, but instead has the form:

$$\sum_{k,l} a^{ij}_{kl} U_{kl} = a^{ij}_{i-1,j}U_{i-1,j} + a^{ij}_{i+1,j}U_{i+1,j} -a^{ij}_{i,j}U_{i,j} + a^{ij}_{i,j-1}U_{i,j-1} + a^{ij}_{i,j+1}U_{i,j+1} = f_{i,j}$$

with non-negative off-diagonal ($a^{ij}_{k,l} \ge 0$ for all off-diagonal entries $(i,j) \ne (k,l)$) and zero row sums, $$\sum_{kl} a^{ij}_{kl} = 0$$ for all $i,j$. From these conditions, the matrix $A = (a^{ij}_{kl})$ is the negative of an $M$-matrix.

The problem is sparse and diagonally dominant (I would characterize it as "weakly diagonally dominant" because the diagonal term is exactly equal to the sum of the off diagonal elements of the same row). Since it exhibits a similar structure to the five-point Laplacian (except perhaps the possibility of being symmetric?), I figured that there may be some similarities between methods to solve this system of equations and those of the five-point discrete Laplacian. Specifically, I was thinking that one of the following may be true:

1. We can compact the sparsity of the matrix into a thinly banded system (since there are only five variables per equation) and use a direct solver.
2. We can use some sort of Fast Fourier Transform to solve it. (I think this one is highly unlikely though.)
3. We can use some sort of multigrid approach to solve the system iteratively.

In the end, the systems will be very large and a parallel implementation will be necessary. Any feedback on fast methods to approach this problem in parallel would be greatly appreciated.

• @Paul Your stencil notation does need correcting. In the general case, you need four indices, in which case the matrix equation reads $\sum_{kl} a_{ij}^{kl} U_{kl} = f_{ij}$ for each grid point $(i,j)$. Then you can ask for zero row sums $\sum_{kl} a_{ij}^{kl} = 0$ and non-negative off-diagonal entries $a_{ij}^{kl} \ge 0$ for $(k,l) \ne (i,j)$. – Jed Brown Jan 19 '12 at 21:32