% Jacobi iteration for solving Ax = b
n = 5;
A = rand(n) + 2*eye(n);
b = rand(n,1);
M = diag(diag(A)); % preconditioner
invM = diag(A).^-1; % vector of diagonal entries of inverse M (instead of inv(M))
% Initial guess
xk = zeros(n,1); % trivial guess, will be updated
rk = b; % first residual is rk = b-A*xk with xk zero vector
resNorm = norm(rk); % norm of the residual to measure convergence gap
resNormVec = []; % vector of norms of rk
maxIter = 1000;
% Iterate:
% Using while loop
while resNorm > 10^-6 & length(resNormVec)