Реализация линейного решателя

Question

Я натолкнулся на следующую реализацию линейного решателя, которая используется для реализации метода стабильных жидкостей Джоса Стама:

/**
 * 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);
    }
}

(Реализация с http://www.multires.caltech.edu/teaching/demos/java/FluidSolver.java)

Я смог собрать x а также x0 являются матрицами, где I это просто простая функция индексации в массиве, и b это просто флаг, для которого устанавливаются граничные условия.

Но я не уверен, что именно делает этот код. Я думал, что линейные решатели в общем пытаются вычислить x за Ax = b, Но в этом коде, похоже, нет векторов, только матрицы.

У кого-нибудь есть идеи, как работает этот код? Я пытаюсь сделать его более эффективным с помощью метода предобусловленных сопряженных градиентов. Но я не понимаю эталонную реализацию здесь.