%function made for the sake of testing

function test_gaussian()
%just input the values here

    A = [ 20, -10,   0; -10,  30, -10; 0, -10,  20 ]; 
    b = [12; 0; -6];
    
    % making use of the function
    Current = gaussianElimination(A, b);
   
    
    
    % just output 
    disp('Final Solution (Currents):');
    disp(Current);
end


function x = gaussianElimination(A, b)
    n = length(b);
    
    % combining A and b into one matrix
    Ab = [A, b];

    % forward elimination & pivoting

    k = 1; % Loop counter 
    while k <= (n - 1)
        
        % partial pivoting
        largestVal = abs(Ab(k, k));
        maxRow = k;
        
        r = k + 1;
        while r <= n
            if abs(Ab(r, k)) > largestVal
                largestVal = abs(Ab(r, k));
                maxRow = r;
            end
            r = r + 1;
        end
        
        % Swap the current row (k) with the row containing the largest value
        if maxRow ~= k %this is not equal
            j = 1;
                while j <= (n + 1)
                tempVal = Ab(k, j);
                Ab(k, j) = Ab(maxRow, j);
                Ab(maxRow, j) = tempVal;
                j = j + 1;
                end
                end
        
        pivot = Ab(k, k);
        if abs(pivot) < 1e-12
            error('Error: Division by zero or singular matrix detected.');
        end
        % Eliminating rows
        i = k + 1;
        while i <= n
            factor = Ab(i, k) / pivot;
            
            % Subtract pivot row from the current row
            j = k;
            while j <= (n + 1)
                Ab(i, j) = Ab(i, j) - (factor * Ab(k, j));
                j = j + 1;
            end
            i = i + 1;
        end
        k = k + 1;
    end
    % back substiution
    % creating empty vector to store (currents)
    x = zeros(n, 1);
    i = n;
    while i >= 1
        runningSum = 0;
        j = i + 1;
        while j <= n
            runningSum = runningSum + (Ab(i, j) * x(j));
            j = j + 1;
        end
        x(i) = (Ab(i, end) - runningSum) / Ab(i, i); %end has special function of getting the last index of a function
        i = i - 1;
    end
end