Gauss

                Never    
C++
       
#include<bits/stdc++.h>
using namespace std;

vector<vector<double>> row_echelon(vector<vector<double>> A)
{
    int n = A.size();
    int m = A[0].size();
    int i = 0, j = 0;
    while(i < n && j < m)
    {
        int k = i;
        while(k < n && A[k][j] == 0)
            k++;
        if(k == n)
        {
            j++;
            continue;
        }
        swap(A[i], A[k]);
        double d = A[i][j];
        for(int l = 0; l < m; l++)
            A[i][l] /= d;
        for(int k = i + 1; k < n; k++)
        {
            double d = A[k][j];
            for(int l = 0; l < m; l++)
                A[k][l] -= d * A[i][l];
        }
        i++;
        j++;
    }
    return A;
}

vector<vector<double>> reducedRowEchelon(vector<vector<double>> matrix) 
{
    int row=matrix.size();
    int col=matrix[0].size();
    int lead=0;
    while(lead<row)
    {
        double d,m;
        for(int r=0;r<row;r++)
        {
            d=matrix[lead][lead];
            m=matrix[r][lead]/matrix[lead][lead];
            for(int c=0;c<col;c++)
            {
                if(r==lead)
                    matrix[r][c]/=d;
                else
                    matrix[r][c]-=matrix[lead][c]*m;
            }
        }
        lead++;
    }
    return matrix;
}


vector<double> gauss_jordan(vector<vector<double>> A,vector<double> B)
{
    vector<vector<double>> C = row_echelon(A);
    int n = C.size();
    int m = C[0].size();
    vector<double> X(n);
    for(int i = n - 1; i >= 0; i--)
    {
        double sum = 0;
        for(int j = i + 1; j < n; j++)
            sum += C[i][j] * X[j];
        X[i] = C[i][n] - sum;
    }
    return X;
}

int main()
{
    int n;
    cin >> n ;
    vector<vector<double>> A(n, vector<double>(n+1));
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
            cin >> A[i][j];
    vector<double> B(n) ;
    
    for(int i = 0; i < n; i++)
    {
        cin >> A[i][n];
    }
    //vector<vector<double>> C = reducedRowEchelon(A);
    vector<vector<double>> C = row_echelon(A);
    cout<<"Row echelon form:"<<endl;
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n+1; j++)
            cout << C[i][j] << " ";
        cout << endl;
    }
    cout << endl;
    vector<double> X = gauss_jordan(A,B);
    cout<<"Solution:"<<endl;
    for(int i = 0; i < n; i++)
        cout << X[i] << " ";
    cout << endl;
    
    return 0;
}
// 3
// 1 2 3
// 4 5 6
// 7 8 9

vector<vector<double>> inverse(vector<vector<double>> A)
{
    int n = A.size();
    vector<vector<double>> B(n, vector<double>(2 * n));
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < 2 * n; j++)
        {
            if(j < n)
                B[i][j] = A[i][j];
            else
                B[i][j] = (i == j - n);
        }
    }
    vector<vector<double>> C = row_echelon(B);
    vector<vector<double>> D(n, vector<double>(n));
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
            D[i][j] = C[i][j + n];
    return D;
}

Raw Text