%
% gaussbs.m - a function that computes the solutions of a
%             system of linear equations Ax=b  for given
%             a nxn matrix A and a vector b if the solution
%             is unique.  If the solution is not uniquet, 
%             a message will be sent out.
%             
% To call it:  x=gaussbs(A,b,n) or
%              [x,umat]=gaussbs(A,b,n)
%
%         n - size of the problem
%         x - an nx1 vector contains solutions.
%         A - an nxn coefficient matrix.
%         b - an nx1 right-hand side vector.
%         umat - an upper triangular matrix 
%                where Ux=v is equivalent to Ax=b
%         ct - a positive integer 
%              which indicates the number of type I elementary row operations
%              used during Gaussian elimination
%  
% by M. Chen
function [xv,umat,ct]=gaussbs(A,b,n);
amat=[A b];
ct=0;
%
% Begin the Gaussian Elimination
%
flag=1;
for i=1:n-1,
    if amat(i,i)==0,
       flag=0;
       k=i+1;
       while flag==0,
             if amat(k,i)~=0,
                flag=1;
                tempv(i:n+1)=amat(i,i:n+1);
                amat(i,i:n+1)=amat(k,i:n+1);
                amat(k,i:n+1)=tempv(i:n+1);
                ct=ct+1;
             else
                k=k+1;
                if k>n,
                   disp('The system has no unique solution!')
                   return
                end;
             end;
       end;
    end;
    if flag==1,
       for j=i+1:n,
           m=-amat(j,i)/amat(i,i);
           amat(j,i+1:n+1)=amat(j,i+1:n+1)+m*amat(i,i+1:n+1);
       end;
    end;
end;
%
% the end of Gaussian Elimination and 
% the beginning of the Backward Substitutuion
%
if flag==1
   if amat(n,n)==0,
      disp('The system has no unique solution! ')
      flag=0;
      return
   end;
   if flag==1,
      xv(n,1)=amat(n,n+1)/amat(n,n);
      for i=n-1:-1:1,
          ss=amat(i,n+1);
          for j=i+1:n,
              ss=ss-amat(i,j)*xv(j);
          end;
          xv(i,1)=ss/amat(i,i);
      end;
      disp('The solution vector is: ')
      umat=zeros(n);
      for ii=1:n;
          umat(ii,ii:n)=amat(ii,ii:n);
      end;
   end
end;