clear; clc;

n=5;
a=rand(5,5); aa=a; %   a=triu(aa);  aa=a
b=rand(5,1);  bb=b;

for ic=1:(n-1)   %..... loop over all columns except last
    for ir=(ic+1):n        %...... loop over rows below the diagonal
        mult=a(ir,ic)/a(ic,ic);
        for icc=ic:n
            a(ir,icc) = a(ir,icc) - mult*a(ic,icc);  %  subtract multiple of above row from this row
        end
        b(ir) = b(ir) - mult*b(ic);  %  do the same with the right hand side
    end
end
% a, return;  %  a should be all zeros below the diagonal now

 %.... back substitution
x(n)=b(n)/a(n,n);
for ir=(n-1):-1:1  %   backwards over rows
    for ic=(ir+1):n
        b(ir) = b(ir) - a(ir,ic)*x(ic); % subtract the known part
    end   
    x(ir)=b(ir)/a(ir,ir);  %   and divide by the diagonal
end

for ir=1:n;   %  multiply back with the original matrix
    b2(ir)=0;
    for ic=1:n
        b2(ir) = b2(ir) + aa(ir,ic)*x(ic);
    end
end
b3=aa*x';    %  backwards in matlab matrix notation

fprintf('bb: ');fprintf(' %20.13f ',bb);  fprintf('\n');
fprintf('b2: ');fprintf(' %20.13f ',b2);  fprintf('\n');
fprintf('b3: ');fprintf(' %20.13f ',b3);  fprintf('\n\n');

x2=aa\bb;
fprintf('x:  ');fprintf(' %20.13f ',x );  fprintf('\n');
fprintf('x2: ');fprintf(' %20.13f ',x2);  fprintf('\n');


