%Test script.
%Homework 9.
%MA 321 Fall 2004.
%Russell Brown


trials = 10;
nrows = 40; 

error = zeros(trials,3);
normResidual = zeros(trials,3);
one = ones(nrows,1);

for  j=1:trials
 
 a = rand(nrows,nrows);;

 b = a * one;
 x = mygauss(a,b);
 error(j,1) = norm(x-one);
 residual = a*x- a*one;
 normResidual(j,1) = norm(residual);

 [c index] = mygauss2(a);
 x = mysolve(c, index, b);
 
 error(j,2) = norm(x-one);
 residual = a*x -a*one;
 normResidual(j,2) = norm(residual);
 
 b = a * one;
 x = a\b;
 
 error(j,3) = norm (x-one);
 residual = a*x - b ;
 normResidual(j,3) = norm(residual);
end

disp('Error for random matrices.') 
disp('Naive Gauss, Gauss with pivoting, Matlab.')
format short e
error/norm(one)
normResidual


%Test on vandermonde matrices.

trials = 14;
error = zeros(trials,3);

for i = 1:trials
    
    vec = [1:1:i];
    a = vander(vec);
    one = ones(i,1);
    b = a*one;
    
    x = mygauss(a,b);
    error(i,1) = norm(x-one);
    normResidual(i,1) = norm(a*x -b);
    
    [c index] = mygauss2(a);
    x = mysolve(c, index, b);
    
    error(i,2) = norm(x-one);
    normResidual(i,2) = norm(a*x-b);
    
    x = a\b;
    
    error(i,3) = norm(x-one);
    normResidual(i,3) = norm(a*x-b);
    
    error(i,[1:3]) = error(i,[1:3])/norm(one);
    
end;
    
    
disp('Error for van der Monde matrices.')
disp('Naive Gauss, Gauss with pivoting, Matlab.')

error
normResidual

%Test which is cooked to illustrate the advantage of pivoting. 

trials = 10;
nrows = 20; 

error = zeros(trials,3);
normResidual = zeros(trials,3);
one = ones(nrows,1);

for  j=1:trials
 %Set up a. 
 a = eye(nrows);
 a = a([nrows:-1:1], [1:nrows]);
 a = a + 10^(-j)*rand(nrows,nrows);;

 b = a * one;
 x = mygauss(a,b);
 error(j,1) = norm(x-one);
 residual = a*x- a*one;
 normResidual(j,1) = norm(residual);


 [c index] = mygauss2(a);
 x = mysolve(c, index, b);
 %Debugging suggestion: examine the index vector. It should be [nrows:-1:1]
 
 error(j,2) = norm(x-one);
 residual = a*x -a*one;
 normResidual(j,2) = norm(residual);
 
 b = a * one;
 x = a\b;
 
 error(j,3) = norm (x-one);
 residual = a*x - b ;
 normResidual(j,3) = norm(residual);
end

disp('Error for cooked examples.') 
disp('Naive Gauss, Gauss with pivoting, Matlab.')

format short e
error/norm(one)
normResidual
   