% Matlab code testgecp.m
% For "Applied Numerical Linear Algebra",  Question 2.14
% Written by James Demmel, Sept 25 2004
%
% Test Gaussian elimination with complete pivoting
%
testcase = 0; clear results;
% reset random number generator to initial state
% this ensures that same matrices are tested each time.
randn('state',0);
%
% generate random matrices with condition number 10^logcond
for logcond = 1:3:13;
   % generate random matrices with dimensions n
   for n = 20:40:100
      % generate one random matrix of each condition number and dimension
      A = randmatgen(logcond,n);
      % generate random solution x to A*x=b
      x = randn(n,1);
      b = A*x;
      % test it
      testcase = testcase + 1;
      results(testcase,:) = testge(A,x,b);
   end
end
%
% generate special matrices with large pivot growth
for n = 30:10:80
   A = diag(ones(n,1)+rand(n,1)*1e-6) - tril(ones(n,n),-1);
   A(:,n) = ones(n,1) + rand(n,1)*1e-6;
   % generate random solution x = A*x=b
   x = randn(n,1);
   b = A*x;
   % test it
   testcase = testcase + 1;
   results(testcase,:) = testge(A,x,b);
end
format short e
disp('See documention of testge.m for more complete description of results')
disp(' ')
disp('Results(1) = dimension of A')
disp('Results(2) = log10(cond(A))')
disp('     Testcase   Results(1)   Results(2)')
[(1:testcase)',results(:,1:2)]
disp(' ')
disp('Results(3) = backward error of factorization from gepp')
disp('Results(4) = pivot growth from gepp')
disp('Results(5) = backward error of factorization from gecp')
disp('Results(6) = pivot growth from gecp')
disp('Results(7) = backward error of solution from gepp')
disp('Results(8) = backward error of solution from gecp')
disp('     Testcase   Results(3)   Results(4)   Results(5)   Results(6)   Results(7)   Results(8)')
[(1:testcase)',results(:,3:8)]
disp(' ')
disp('Results(9)  = condition estimate from gepp')
disp('Results(10) = condition estimate from gecp')
disp('Results(11) = explicit computation of condition number')
disp('     Testcase   Results(9)  Results(10)  Results(11)')
[(1:testcase)',results(:,9:11)]
disp(' ')
disp('Results(12) = true error in x computed by  gepp')
disp('Results(13) = error bound on x computed by gepp')
disp('Results(14) = true error in x computed by gecp')
disp('Results(15) = error bound on x computed by gecp')
disp('     Testcase  Results(12)  Results(13)  Results(14)  Results(15)')
[(1:testcase)',results(:,12:15)]