%% Function iterSolve
% Solves a matrix system Ax=b using either Jacobian iterations or Gauss-
% Seidel interations. Matrix must be inverible and the 1 Norm of the matrix 
% must be less than 1, otherwise code will not run.
%
%EXAMPLES OF USE
%   iterSolve(A,b,opt);
%   [x,count] = iterSolve(A,b,opt);
%   [x,count,dxVec] = iterSolve(A,b,opt);
%   [x,count,dxVec,err] = iterSolve(A,b,opt);
%
%TEST EXAMPLE CODE
% opt.maxIt = 20000;
% opt.errTol = 1.e-8
% for k = 2:7
%     n = 2^k;
%     fprintf('Test %1.0f: A is an %3.0f by %3.0f matrix \n',k-1,n,n)
%     D=sparse(1:n,1:n, 3*ones(1,n),n,n);
%     E=sparse(2:n,1:n-1, -ones(1,n-1),n,n);
%     A=1/n*(D + E + E');%create text matrix(tridiagonal)
%     b=rand(n,1);
%     opt.method = 1;%Jacobi
%     [x,count,dxVec,err] = iterSolve(A,b,opt);
%     figure(k-1)
%     loglog(1:count,err,'LineWidth',2)
%     grid on
%     hold on
%     opt.method = 2;%GS
%     [x,count,dxVec,err] = iterSolve(A,b,opt);
%     loglog(1:count,err, '*r')
%     grid on
%     title('Jacobi Iterations: Log plot Count VS max|A^{-1}b - x^k| ');
%     xlabel('Count');
%     ylabel('Error');
%     legend('Jacobi Iterations','Gauss-Seidel Iterations')
%     hold off
%     fprintf('Number of iterations: %1.0f: with err: %1.4e \n',count,err(end))
% end
%
%INPUT
%   A = Square matrix
%   b = solutions vector
%   opt is an optional data structure with the following elements
%       opt.guess = initial guess (zero vector = default)
%       opt.errTol = error tolerance (1.e-10 = default)
%       opt.maxIt = max number of iterations run (1000 = default)
%       opt.method = 1 Jacobi iterations
%                    2 Gauss Seidel iterations (default)
%   
%OUTPUT
%   x = solutions of the system Ax=b
%   count = number of iterations
%   dxVec = vector of residuals for each iteration
%   err = vector of max(abs(sol - x)) for each iteration

function [x,count,dxVec,normDxVec] = iterSolve(A,b,opt)
%
n = length(b);
%
if nargin == 2
   opt = [];
end
%
if isfield(opt,'errTol') == 0
    errTol = 1.e-10;
else
    errTol = opt.errTol;
end
%
if isfield(opt,'guess') == 0
    x = zeros(n,1);
else
    x = opt.guess;
end
%
if isfield(opt,'method') == 0
    method = 2;
else
    method = opt.method;
end
%
if isfield(opt,'maxIt') == 0
    maxIt = 1000;
else
    maxIt = opt.maxIt;
end
if norm(A,1) > 1 % Check if matrix is "good"
    disp('1-norm of matrix is greater than 1, recondition matrix');
    x = [];
    count = [];
    dxVec = [];
    normDxVec = [];
    return;
elseif det(A) == 0
    disp('Matrix is near singular');
    x = [];
    count = [];
    dxVec = [];
    normDxVec = [];
    return;
end
dxVec = zeros(n,maxIt);
err = zeros(1,maxIt);
count = 0;
sol = A\b;
switch method
    case 1
        while max(abs(sol - x))  > errTol && count < maxIt %Jacobi iterations
            R = b - A*x; %construct residual
            dx = R./diag(A); %find step size
            x = x + dx; %update iteration
            count = count + 1; %iteration number
            dxVec(:,count) = dx; %record residuals
            err(count) = max(abs(sol - x));%record errors
        end
    case 2
        while max(abs(sol - x))  > errTol && count < maxIt %GS iterations
            R = b - A*x; %construct residual
            dx = tril(A)\R; %find step size
            x = x + dx; %update iteration
            count = count + 1; %iteration number
            dxVec(:,count) = dx; %record residuals
            err(count) = max(abs(sol - x));%record errors
        end
end
dxVec = dxVec(:,1:count);
normDxVec = err(1:count);