Решение СЛАУ методом простых итераций

Нужно сделать программу, которая будет решать СЛАУ методом простых итераций. Есть две вектора, один - матрица с корнями, второй - со свободными коэффициентами. Когда программа считает первую итерацию - выводит правильный ответ и корни хранит в векторе, хочу чтобы в последующих итерациях она брала корень из этого вектора и вычисляла, но так же нужно, чтобы эти корни не повторялись. т.е. если считать итерацию для второго корня, нужно чтобы она брала 0,1 и 3 корни, не могу понять как реализовать цикл.

for (int i = 0; i < matrix.size(); i++)
            {
                if (j > 0)
                {
                    X[i] = free[i] / matrix[i][i];
                    if (g > 0)
                    {
                        for (int j = 1; j < matrix.size(); j++)
                        {
                            if (i == j)
                            {
                                continue;
                            }
                            temp[i] += matrix[i][j-i]*X[j-i];
                        }
                        X[i] = (free[i] - temp[i]) / matrix[i][i];
                    }
                    cout << fixed << setprecision(3) << X[i] << " ";
                }
            }

Вот весь код:

int main()
{
    float eps = 1000;
    setlocale(LC_ALL, "Russian");
    float M=0.89, N=0.08, P=-1.21;
    vector<vector<float>>matrix = { {M,-0.04,0.21,-18},
                                 {0.25,-1.23,N,-0.09},
                                 {-0.21,N,0.8,-0.13},
                                 {0.15,-1.31,0.06,P} };
    vector<float> free = { -1.24,P,2.56,M };
    cout << "Исходная матрица:" << endl;
    for (int i = 0; i < matrix.size(); i++)
    {
        for (int j = 0; j < matrix.size(); j++)
        {
            cout << matrix[i][j] << "\t";
        }
        cout << free[i] << endl;;
    }
    vector<float> X = { 0,0,0,0};
    vector<float> temp = { 0,0,0,0 };
    cout << "Введите количество итераций: ";
    int n;
    cin >> n;
    int j = 4;
    for (int g = 0; g < n; g++)
    {
        cout << g + 1 << "\t";
        for (int i = 0; i < matrix.size(); i++)
        {
            if (j > 0)
            {
                X[i] = free[i] / matrix[i][i];
                if (g > 0)
                {
                    for (int j = 1; j < matrix.size(); j++)
                    {
                        if (i == j)
                        {
                            continue;
                        }
                        temp[i] += matrix[i][j-i]*X[j-i];
                    }
                    X[i] = (free[i] - temp[i]) / matrix[i][i];
                }
                cout << fixed << setprecision(3) << X[i] << " ";
            }
        }
        cout << endl;
    }

}