Linear solver implementation

Question

I came across the following implementation of a linear solver, which is used to implement Jos Stam's stable fluids method:

/**
 * Iterative linear system solver using the Gauss-sidel
 * relaxation technique. Room for much improvement here...
 **/
void linearSolver(int b, float[] x, float[] x0, float a, float c)
{
    for (int k = 0; k < 20; k++)
    {
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= n; j++)
            {
                x[I(i, j)] = (a * ( x[I(i-1, j)] + x[I(i+1, j)]
                                +   x[I(i, j-1)] + x[I(i, j+1)])
                                +  x0[I(i, j)]) / c;
            }
        }
        setBoundry(b, x);
    }
}

(Implementation from http://www.multires.caltech.edu/teaching/demos/java/FluidSolver.java)

I have been able to gather x and x0 are matrices, where I is just a simple indexing function into the array, and b is simply a flag for which boundary conditions to set.

But I'm not sure what exactly this code is doing. I thought linear solvers in general try to compute x for Ax = b. But there don't seem to be any vectors in this code, only matrices.

Does anyone have any idea how this code works? I am attempting to make it more efficient via the method of preconditioned conjugate gradients. But I don't understand the reference implementation here.