function iter = myjacobi (myn,nproc,tol,maxiter)
%%%%%%%%%%%%%%%%%%%
%%% 
%%%%%%%%%%%%%%%%%%%%
%%% calling sequence  .... myjacobi(20,1,1e-3,1000)
%%% myjacobi(10,2,1e-3,1000)

%%%%% PAR 
%% every processor needs to know his number <p> as well nproc, and myn,
%% only proc 1 needs to know tol and maxiter 
%%%%%

n = myn*nproc; %% size of global problem

%%% 
h = 1/(n+1);    %% n,h are global 

%%% DATA STRUCTURES
%%% we need b, A, unew, Aleft, Aright, uleft, uright
myx = zeros(myn,nproc);
for p=1:nproc
    myx(1:myn,p) = (p-1)*myn*h+(1:myn)'*h; 
    myb = h * h * ones(myn,1);
end;

%%% cheating ... matrix is the same on all processors 
myA = zeros (myn,myn);
for i = 1:myn myA(i,i)=2;end;
for i = 2:myn myA(i,i-1)=-1;end;
for i = 1:myn-1 myA(i,i+1)=-1;end;
%%% off-block diagonal matrix proc components   
Aright = -1;
Aleft = -1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% this is for illustraiton
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% only
%%% global matrix: for illustration only 
A = zeros(n,n);
for i = 1:n A(i,i)=2;end;
for i = 2:n A(i,i-1)=-1;end;
for i = 1:n-1 A(i,i+1)=-1;end;
x = (1:n)*h;    
b = h * h * ones(n,1);
utrue = A \ b;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% read or communicate initial guess, perhaps different on every processor
%% cheating here ... it will be the same on every processor
u0 = zeros(myn,1);

%% initialize guess on every processor
for p =1:nproc u(:,p) = u0;end;

for iter = 1:maxiter  

    for p = 1:nproc
        
        %% every processor has to have current values in ghost nodes 
        %%% SEND    your values
        %%% RECEIVE your values 
        %%% here this is done by grabbing elements from global memory
        
        if p > 1        uleft(p) = u(myn,p-1);end;
        if p < nproc    uright(p) = u(1,p+1);end;
        
    end;
        
    for p = 1:nproc
        %% compute current residual and its norm  
        myr = myb - myA*u(:,p);          

        %% do this in parallel: use ghost cells
        if p > 1 myr(1) = myr(1) - Aleft*uleft(p);end;
        if p < nproc myr(myn) = myr(myn) - Aright*uright(p);end;
  
        rerr(p) = norm(myr);         %% local calculation 
    end;
        
    
    %% collect all values of local calculations 
    s = 0; for p = 1:nproc s = s + rerr(p)*rerr(p);end;
    glob_err = sqrt(s);    
    
    
    %%% only processor 1 communicates with the world 
    if p == 1 
        fprintf('%d %g\n',iter,glob_err); 
    end;
    
    %%% every processor must know that it is time to stop 
    if glob_err < tol break;end;    
   
    %%%% perform actual loop on every processor

    for p = 1:nproc
        
        for k = 1:myn
            %% collect right hand side terms for every row 
            s = 0;
            for j = 1:k-1 s = s + myA(k,j)*u(j,p);end;
            for j = k+1:myn s = s + myA(k,j)*u(j,p);end;
        
            %%% incorporate values from left and right 
            if k ==1 if p > 1       s = s + Aleft*uleft(p);end;end;
            if k == myn if p < nproc   s = s + Aright*uright(p); end; end;

            %% compute new guess in row k
            unew (k,p) = (myb(k) - s)/myA(k,k);
        end;
        
        u(:,p) = unew(:,p);
    end
        
    %% visualization of current guess and true solution 
    plot(x,utrue);hold on;for p=1:nproc plot(myx(:,p),unew(:,p)); axis([0 1 0 1]);   end; hold off;pause(0.1);
        
end

    
