function [x,resvec]=jacobi(A,b,x,omega,n)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% relaxed Jacobi iteration to solve Ax=b:
%
%         D (xnew-x) = omega (b-Ax)   where  D= diag(A)
%
%  [x,resvec]=jacobi(A,b,x,omega,n) specifies initial guess x, omega
%  and number of iterations n. Returns "x" and "resvec" the relative
%  residual history.
%
%  [x,resvec]=jacobi(A,b)   picks x=0, omega=2/3, n=30
%
%  [x,resvec]=jacobi(A,b,x) uses "x" as guess, omega=2/3, n=30
%
%
%  Fabian Waleffe, 11/15/2000 demo for Math 801
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize
  
if nargin<3 
   x=zeros(size(b)); omega=2/3; n=30;
end 
if nargin<4 
   omega=2/3; n=30;
end 

  D=diag(A);
  r= b-A*x; rTr = r'*r; bTb = b'*b;   

  relres= sqrt(rTr/bTb);
  resvec= [relres];

%% Jacobi iterations

 for i=1:n,
    
    x = x + omega*r./D;     
    r = b - A*x;       

    relres=  sqrt(r'*r/bTb);
    resvec=[resvec;relres];
  end



