%risoluzione di un sistema triangolare -con pivoting
clear all
clc

A=[0.0001 0.5;0.4 0.3];
b=[0.5;0.1];
[A,b] = GAUSSMIO(A,b);
[n,~]=size(A);

%A questo punto posso applicare risoluzione per sistema triangolare
%superiore
x=zeros(n,1);
x(n)=b(n)/A(n,n);
s=0;
for i=n-1:-1:1
    for j=i+1:n
        s=s+A(i,j)*x(j);
    end
    x(i)=(b(i)-s)/A(i,i);
end

A=[0.0001 0.5;0.4 0.3];
b=[0.5;0.1];

[A,b] = GaussJordan(A,b);
%soluzione di un sistema diagonale

for i=1:n
    x(i)=b(i)/A(i,i);
end
x

function [A,b] = GAUSSMIO(A,b)
%metodo di Gauss
[n,~]=size(A);

for k=1:n-1
    [~,r]=max(A(k:n,k))
    temp=A(k,:);
    A(k,:)=A(r,:);
    A(r,:)=temp;

    temp=b(k);
    b(k)=b(r);
    b(r)=temp;

    for i=k+1:n
        mik=-A(i,k)/A(k,k);
        b(i)=b(i)+mik*b(k);

        for j=k+1:n
            A(i,j)=A(i,j)+mik*A(k,j);
        end
        
    end
end
A=triu(A);
end


function [A,b] = GaussJordan(A,b)

%metodo di Gauss
[n,~]=size(A);

for k=1:n-1
    [~,r]=max(A(k:n,k));
    temp=A(k,:);
    A(k,:)=A(r,:);
    A(r,:)=temp;

    temp=b(k);
    b(k)=b(r);
    b(r)=temp;

    for i=k+1:n
        mik=-A(i,k)/A(k,k);
        b(i)=b(i)+mik*b(k);

        for j=k+1:n
            A(i,j)=A(i,j)+mik*A(k,j);
        end
        
    end
end

A=triu(A);

for k=n:-1:2
    [~,r]=max(A(i:k,k));
    temp=A(k,:);
    A(k,:)=A(r,:);
    A(r,:)=temp;

    temp=b(k);
    b(k)=b(r);
    b(r)=temp;
    for i=k-1:-1:1
        mik=-A(i,k)/A(k,k);
        b(i)=b(i)+mik*b(k);

        for j=n:-1:k
            A(i,j)=A(i,j)+mik*A(k,j);
        end
    end
end

end