% Compile test results for Programming Assignment 3, Ma128a
%
% This program solves randomly generated linear systems using
% three algorithms:
%   GENP = Gaussian elimination with no pivoting
%   GEPP = Gaussian elimination with partial pivoting
%   GECP = Gaussian elimination with complete pivoting
%
% For each test matrix, this program
%   1) Solves Ax=b using the 3 algorithms GENP, GEPP, and GECP.
%      Some of the test matrixes are very hard, and Matlab
%      will print a lot of warning messages that you should ignore.
%      It also gets the pivot growth PG, defined as
%      the largest entry in L times the largest entry in U divided by
%      the largest entry in A.
%      PG should be at least about 1, and we expect an inaccurate answer if PG>>1.
%
%   2) Computes the error 
%          Err = norm(x_true - x_computed)/norm(x_true)
%      for each algorithm. 
%      When Err is around machine epsilon, 1e-16,
%      this means x_computed is about as accurate as possible
%      on the machine.
%      When Err is around 1 or greater, this means 
%      x_computed is completely inaccurate.
%
%   3) Computes an error bound Ebnd, which should bound Err
%      if the algorithm does not experience "pivot growth", i.e. PG>>1.
%      Otherwise, PG*Ebnd should bound Err.
%
%   4) Computes the ratio R = Err/Ebnd for each algorithm.
%      R should be less than 1 (or perhaps a little larger)
%      as long as there is little pivot growth, i.e. when
%      the implementation of Gaussian elimination is working
%      as well as can be expected.
%      Also, we expect R to be not much smaller than 1,
%      say .01 or greater, If it is smaller, the algorithm was 
%      more accurate than expected, a lucky accident of rounding error.
%
%   5) Compute the ratio R_PG = Err/(Ebnd*PG) for each algorithm.
%      This should be less than 1 (or perhaps a little larger) even
%      if there is pivot growth.
%
%   6) Prints a table of the above results
%
disp('tstGE - test driver for Program 3, Math 128a')
disp('Ignore messages: "Warning: Matrix is close to singular or badly scaled."')
clear result
i=0;
N = [5,10,25,50,100];    %  dimensions to test
lengthN = length(N);     %  number of dimensions to test
numcases = 5;
for ncase = 1:numcases    %  for all 5 kinds of test matrices
   disp(['Starting case = ',int2str(ncase),' ...'])
   for n = N                   %  for the matrix dimensions in N
     i = i+1;                         % increment example #
     [A,x,b]= matgen(ncase,n);        % generate a test problem A*x=b
     [x_NP,PG_NP] = genpsolve(A,b);   % solve using GENP, get pivot growth
     [x_PP,PG_PP] = geppsolve(A,b);   % solve using GEPP, get pivot growth
     [x_CP,PG_CP] = gecpsolve(A,b);   % solve using GECP, get pivot growth
     Err_NP = norm(x_NP-x)/norm(x);   % compute GENP error
     Err_PP = norm(x_PP-x)/norm(x);   % compute GEPP error
     Err_CP = norm(x_CP-x)/norm(x);   % compute GECP error
     Ebnd = eps*cond(A);              % compute error bound
     R_NP = Err_NP/Ebnd;              % compute error bound ratio for GENP
     R_PP = Err_PP/Ebnd;              % compute error bound ratio for GEPP
     R_CP = Err_CP/Ebnd;              % compute error bound ratio for GECP
     R_PG_NP = Err_NP/Ebnd/PG_NP;     % compute error bound ratio for GENP,
                                      %  including pivot growth
     R_PG_PP = Err_PP/Ebnd/PG_PP;     % compute error bound ratio for GEPP
                                      %  including pivot growth
     R_PG_CP = Err_CP/Ebnd/PG_CP;     % compute error bound ratio for GECP
                                      %  including pivot growth
%    Collect results
     result(1:16,i)= [i,ncase,n,Err_NP,Err_PP,Err_CP,PG_NP,PG_PP,PG_CP, ...
                      Ebnd,R_NP,R_PP,R_CP,R_PG_NP,R_PG_PP,R_PG_CP]';
   end
end
% Display results in a nice format
disp(' ')
disp('Result of GENP, GEPP and GECP on test problems A*x=b')
disp(' ')
if (length(N) >= 6 | length(N) <= 3),
    numperrow = 6;
else
    numperrow = length(N);
end
display_data  = num2str(result(:,:),2);
[rdd,cdd]=size(display_data);
for j=1:ceil(numcases*lengthN/numperrow),
   j1 = (j-1)*9*numperrow+1;  % first case to print on this line
   j2 = min([j*9*numperrow,numcases*9*lengthN-1,cdd]); % last case to print
   disp(['Example #  ',display_data(1,j1:j2)])   % Example # 
   disp(['Case       ',display_data(2,j1:j2)])   % Case (matrix type)
   disp(['Dimension  ',display_data(3,j1:j2)])   % Matrix dimension
   disp(['Err_NP     ',display_data(4,j1:j2)])   % error from GENP
   disp(['Err_PP     ',display_data(5,j1:j2)])   % error from GEPP
   disp(['Err_CP     ',display_data(6,j1:j2)])   % error from GECP
   disp(['PG_NP      ',display_data(7,j1:j2)])   % pivot growth from GENP
   disp(['PG_PP      ',display_data(8,j1:j2)])   % pivot growth from GEPP
   disp(['PG_CP      ',display_data(9,j1:j2)])   % pivot growth from GECP
   disp(['Ebnd       ',display_data(10,j1:j2)])  % error bound
   disp(['R_NP       ',display_data(11,j1:j2)])  % error ratio for GENP
   disp(['R_PP       ',display_data(12,j1:j2)])  % error ratio for GEPP
   disp(['R_CP       ',display_data(13,j1:j2)])  % error ratio for GECP
   disp(['R_PG_NP    ',display_data(14,j1:j2)])  % error ratio for GENP
                                                 %  including pivot growth
   disp(['R_PG_PP    ',display_data(15,j1:j2)])  % error ratio for GEPP
                                                 %  including pivot growth
   disp(['R_PG_CP    ',display_data(16,j1:j2)])  % error ratio for GECP
                                                 %  including pivot growth
   disp(' ')
end