%
% GE.m
%
clc; clear all; close all; format long
nseed = fix(rand*1000000); rand('seed',nseed);
%
tic;  costime = toc;
%
disp('  ')
disp('  The program GE.m starts ...')
disp('  ')
%
disp(' To solve a linear system A x = b, then please give A and b...'); disp('  ')
n = input(' the dimension of A = '); if ((n<=0) | ((n-fix(n))~=0) ); error(' n is not an integer!'); end;
index_sol = 1:n; ans_inf_sol = 0;
    ans_1 =  input(' Are A and b created by (1) keyin or (2) random or (3) read-file...');disp('  ');
    if(ans_1 == 1)
          A = zeros(n,n); b= zeros(n,1); ans_ok = 0;
      while (ans_ok ~= 1)
        for jj = 1:n
            for ii =1:n
                A(ii,jj) = input([' A(',num2str(ii),',',num2str(jj),') = ']);
            end
        end
            for ii =1:n
                b(ii) = input([' b(',num2str(ii),') = ']);
            end
         disp(' [ A b ] = '); disp([ A b ]); 
         ans_ok = input(' the data of A and b is correct! (1) Yes; (2) No...');
      end
    elseif(ans_1 == 3)
        create_A_b
        size_A = size(A); size_b=size(b);
        if ( (size_A(1)~=n) | (size_A(2)~=n) | (size_b(1)~=n) | (size_b(2)~=1))
            disp([' size of A is ',num2str(size(A)),';'])
            disp([' size of b is ',num2str(size(b))])            
            error(' size of A or b is incorrect!');
        end
    else
        A = rand(n); b = rand(n,1);
    end
        M = [ A b ];
%
for ii = 1:n
    if (M(ii,ii)~= 0)
            M(ii,ii:end) = M(ii,ii:end)/M(ii,ii);
        for jj = (ii-1):-1:1
            M(jj,ii:end) = M(jj,ii:end) - (M(jj,ii)/M(ii,ii))*M(ii,ii:end);
        end
        for jj = (ii+1):n
            M(jj,ii:end) = M(jj,ii:end) - (M(jj,ii)/M(ii,ii))*M(ii,ii:end);
        end
    else
            ans_ch_r = 0;
        for jj = (ii+1):n
          if (M(jj,ii)~= 0)
              tmpr = M(ii,ii:end); M(ii,ii:end) = M(jj,ii:end); M(jj,ii:end) = tmpr; 
              ans_ch_r = 1; break; 
          end
        end
           if (ans_ch_r == 1)
                 M(ii,ii:end) = M(ii,ii:end)/M(ii,ii);
             for jj = (ii-1):-1:1
                 M(jj,ii:end) = M(jj,ii:end) - (M(jj,ii)/M(ii,ii))*M(ii,ii:end);
             end
             for jj = (ii+1):n
                 M(jj,ii:end) = M(jj,ii:end) - (M(jj,ii)/M(ii,ii))*M(ii,ii:end);
             end
           else
               
                 ans_ch_c = 0;
             for jj = (ii+1):n
               if (M(ii,jj)~= 0)
                  tmpc = M(ii:end,ii); M(ii:end,ii) = M(ii:end,jj); M(ii:end,jj) = tmpc; 
                  ans_ch_c = 1; tmp = index_sol(ii); index_sol(ii) = jj; index_sol(jj) = tmp; break; 
               end
             end
                if (ans_ch_c == 1)
                   M(ii,ii:end) = M(ii,ii:end)/M(ii,ii);
                  for jj = (ii-1):-1:1
                     M(jj,ii:end) = M(jj,ii:end) - (M(jj,ii)/M(ii,ii))*M(ii,ii:end);
                  end
                  for jj = (ii+1):n
                     M(jj,ii:end) = M(jj,ii:end) - (M(jj,ii)/M(ii,ii))*M(ii,ii:end);
                  end
                else
                    if (abs(M(ii,end)) > eps)
                        disp(' [ A  b ] = '); disp([ A  b ]); 
                        disp(' The linear system has no solution!');
                        disp(' GE [ A  b ] = '); disp(M); disp([' with column pivoting index ',num2str(index_sol)]);
                        return
                    else
                        ans_inf_sol = 1;
                    end                    
                end
           end
                   
    end
end

if (ans_inf_sol == 1)
    disp(' [ A  b ] = '); disp([ A  b ]); 
    disp(' The linear system has infinty many solutions!');
    disp(' cp [ A  b ] = '); disp(M); disp([' with column pivoting index ',num2str(index_sol)]);
else
    num_sol = M(:,end);
    residual = norm(A*num_sol - b,inf);
    disp(' [ A x b ] = '); disp([ A num_sol b ]); 
    disp([' the residual is ',num2str(residual)]);
end
%
disp('  ')
disp('  The program is terminated!')
disp([  '  Total is ',num2str(toc-costime),' seconds elapsed time!'])
disp('  ')
%