function x = Gauss(A,b,pd,plot)
%   function x = Gauss_pivot(A,b)
%   Solution of the system of linear equations Ax = B
%   using Gaussian elimination without pivoting
%   Inputs :
%       A is the n-by-n coefficient matrix
%       b is the n-by-1 right-hand-side vector
%       pd, if exist, is the number of significant digits to be used.
%       plot if nonzero, causes a display of elimination.
%   Outputs :
%       x is the n-by-1 solution vector

if nargin<3, pd = 0; end
if nargin<4, plot = 0; end
if isempty(pd), pd = 0; end;
if pd>0,
    A = round_matrix(A,pd);
    b = round_matrix(b,pd);
end
if plot
    figure; axis off;plot_matrix([A b]);pause
end

[n,n1] = size(A);
for i=1:n-1
    m = -A(i+1:n,i)/A(i,i);
    A(i+1:n,:) = A(i+1:n,:)+m*A(i,:);
    b(i+1:n,:) = b(i+1:n,:)+m*b(i,:);
    if pd>0
        A = round_matrix(A,pd);
        b = round_matrix(b,pd);
    end
    if plot
        clf; axis off;
        plot_matrix([A b]);
        pause
    end
end
% back substitution
x(n,:) = b(n,:)./A(n,n);
for i=n-1:-1:1
    x(i,:) = (b(i,:)-A(i,i+1:n)*x(i+1:n,:))./A(i,i);
end
close

function A=round_matrix(A,pd)

[m,n] = size(A);
for ii=1:m
    for jj=1:n
        if abs(A(ii,jj))>10^(-pd)/2
            s = floor(log10(abs(A(ii,jj))));
            fmat = ['%0.' num2str(pd-1) 'f'];
            As = num2str(A(ii,jj)/(10^s),fmat);
            A(ii,jj) = str2num(As)*(10^s);
        else
            A(ii,jj) = 0;
        end
    end
end

function plot_matrix(A)

if ~ishold, hold, end
[m,n] = size(A);
ax = axis;
xlen = ax(2)-ax(1);
ylen = ax(4)-ax(3);
xx = linspace(ax(1)+0.1*xlen,ax(2)-0.1*xlen,n+1);
yy = linspace(ax(3)+0.1*ylen,ax(4)-0.1*ylen,m+1);
for ii=1:n+1;
    plot([xx(ii) xx(ii)],[yy(1) yy(m+1)]);
end
for ii=1:m+1;
    plot([xx(1) xx(n+1)],[yy(ii) yy(ii)]);
end
xx = (xx(1:end-1)+xx(2:end))/2;
yy = (yy(1:end-1)+yy(2:end))/2;
for ii=1:m;
    for jj=1:n;
        hx = text(xx(jj),yy(m+1-ii),num2str(A(ii,jj),'%0.3f'));
        set(hx,'HorizontalAlignment','center');
        set(hx,'VerticalAlignment','middle');
    end
end
axis(ax);
        