can anyone modify this gauss elim code into LU gauss elimina

can anyone modify this gauss elim code into LU gauss elimination

function [ x ] = GaussElim(a,b)

%this function solves a sytem of linear equations [a][x]=[b] using the

%Gaussian elimination method

% Input variables:

% a The matrix of coefficients

% b The right-hand side column vector of constants

% Output variable:

% x A column vector with the solution

ab = [a,b];

[R,C] = size(ab);

for j = 1:R-1

for i = j + 1:R

ab(i,j:C) = ab(i,j:C) - ab(i,j)/ab(j,j)*ab(j,j:C);

end

end

x=zeros(R,1);

x(R) = ab(R,C)/ab(R,R);

for i = R - 1:-1:1

x(i) = (ab(i,C) - ab(i,i+1:R)*x(i+1:R))/ab(i,i);

end

please use comment and write clearly

Solution

#include <iostream>

#include <cmath>

#include <vector>

using namespace std;

void print(vector< vector<double> > A) {

    int n = A.size();

    for (int i=0; i<n; i++) {

        for (int j=0; j<n+1; j++) {

            cout << A[i][j] << \"\\t\";

            if (j == n-1) {

                cout << \"| \";

            }

        }

        cout << \"\ \";

    }

    cout << endl;

}

vector<double> gauss(vector< vector<double> > A) {

    int n = A.size();

    for (int i=0; i<n; i++) {

        // Search for maximum in this column

        double maxEl = abs(A[i][i]);

        int maxRow = i;

        for (int k=i+1; k<n; k++) {

            if (abs(A[k][i]) > maxEl) {

                maxEl = abs(A[k][i]);

                maxRow = k;

            }

        }

        // Swap maximum row with current row (column by column)

        for (int k=i; k<n+1;k++) {

            double tmp = A[maxRow][k];

            A[maxRow][k] = A[i][k];

            A[i][k] = tmp;

        }

        // Make all rows below this one 0 in current column

        for (int k=i+1; k<n; k++) {

            double c = -A[k][i]/A[i][i];

            for (int j=i; j<n+1; j++) {

                if (i==j) {

                    A[k][j] = 0;

                } else {

                    A[k][j] += c * A[i][j];

                }

            }

        }

    }

    // Solve equation Ax=b for an upper triangular matrix A

  vector<double> x(n);

    for (int i=n-1; i>=0; i--) {

        x[i] = A[i][n]/A[i][i];

        for (int k=i-1;k>=0; k--) {

            A[k][n] -= A[k][i] * x[i];

        }

    }

    return x;

}

int main() {

    int n;

    cin >> n;

    vector<double> line(n+1,0);

    vector< vector<double> > A(n,line);

    // Read input data

    for (int i=0; i<n; i++) {

        for (int j=0; j<n; j++) {

            cin >> A[i][j];

        }

    }

    for (int i=0; i<n; i++) {

        cin >> A[i][n];

    }

    // Print input

    print(A);

    // Calculate solution

    vector<double> x(n);

    x = gauss(A);

    // Print result

    cout << \"Result:\\t\";

    for (int i=0; i<n; i++) {

        cout << x[i] << \" \";

    }

    cout << endl;

}