function [x, resvec] = jacobisolve(A, b, niters, viziters)
% jacobisolve : Solve Ax=b by Jacobi iteration
%
% [x, resvec] = jacobisolve(A, b, niters);
%
% Given matrix A and vector b, 
% this runs "niters" iterations of Jacobi to solve for x in A*x=b.
%
% x = jacobisolve(A, b, niters, viziters); 
%  also visualizes x as a solution to the model problem
%  after every "viziters" iterations.
%
% Optionally, [x, resvec] = ... also returns a vector of 
% the norm of the residual b - A*x at each iteration.
%
% The code below follows the demo in CS 111.
%
% John R. Gilbert    13 Jan 2016

n = length(b);

% Split A into diagonal D plus off-diagonal C
d = diag(A);
D = diag(d)
C = A - D

% Initial guess: x = 0
x = zeros(n,1);

resvec = zeros(1,niters);

% Iterate
for iter = 1 : niters
    x = (b - C*x) ./ d
    pause
    res = norm(b - A*x);
    resvec(iter) = res;
    if nargin > 3 && mod(iter,viziters) == 0
        k = sqrt(n);
        surfc(reshape(x,k,k))
        title(sprintf('Iteration %d:  residual = %g, maximum x = %g',...
            iter,res,max(x)));
        pause(.5);
    end
end;
