function z = jacobi2d ( r )

%solves Mz=r by jacobi preconditioner

    n = length ( r );
    m = sqrt(n);
    h = 1.0/(m+1);
    h2 = 1/h^2;
    z1=zeros(m,m);
    z=zeros(n,1);
    u=zeros(m,m);
    u(:)=r;
  
    for i=1:m
        for j=1:m
            %j=i;

            %get the four diagonal bits
            z1(i,j)=(h2/2.0*4.0*A(i,j,m));%*u(i,j);

            %get the other bits -- be sure to handle the edges correctly
            if j~=1-1
              z1(i,j) = z1(i,j) + (h2/2.0*(A(i,j-1,m)));%*u(i,j-1);
            end

            if i~=1-1
              z1(i,j) = z1(i,j) + (h2/2.0*(A(i-1,j,m)));%*u(i-1,j);
            end

            if j~=m+1
              z1(i,j) = z1(i,j) + (h2/2.0*(A(i,j+1,m)));%*u(i,j+1);
            end

            if i~=m+1
              z1(i,j) = z1(i,j) + (h2/2.0*(A(i+1,j,m)));%*u(i+1,j);
            end

            z1(i,j)  = 1/z1(i,j) * u(i,j);
        end
    end
    
    z=z1(:);