%
% gaussbsmr.m - a function that computes the solutions of a
%             system of linear equations with r right side vectors 
%             Ax=B  where 
%             A is an nxn matrix and B is an nxr matrix.
%             If the solutions do not exist,  a message will be printed out.
%             
% To call it:  X=gaussbs(A,B,n,r) 
%
%         n - size of the problem
%         r - number of right side vectors
%         X - an nxr matrix containing solutions.
%         A - an nxn coefficient matrix.
%         B - an nxr matrix containing r right-hand side vectors.
%  
% by M. Chen
function X=gaussbsmr(A,B,n,r);
amat=[A B];
%
% 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+r)=amat(i,i:n+r);
                amat(i,i:n+r)=amat(k,i:n+r);
                amat(k,i:n+r)=tempv(i:n+r);
             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+r)=amat(j,i+1:n+r)+m*amat(i,i+1:n+r);
       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,
      for k=1:r
          xv(n,k)=amat(n,n+k)/amat(n,n);
          for i=n-1:-1:1,
              ss=amat(i,n+k);
              for j=i+1:n,
                  ss=ss-amat(i,j)*xv(j,k);
              end;
              xv(i,k)=ss/amat(i,i);
          end;
      end;
      disp('The solutions are: '),xv
      X=xv;
   end;
end;