%----------------------------------------------------------------------------
%  consider 
%      4x - y + z = 7   =>  x = (7 + y - z)/4
%      4x -8y + z =-21  =>  y = (-21 + 4x -z)/-8
%     -2x + y +5z = 15  =>  z = (15 +2x - y)/5
%
%      i.e  x2 = X(1) = [B(1) - A(1,2)*y - A(1,3)*z]/A(1,1)
%           y2 = X(2) = [B(2) - A(2,1)*x - A(2,3)*z]/A(2,2)
%           z2 = X(3) = [B(3) - A(3,1)*y - A(3,2)*y]/A(3,3)
%      in general, 
%         let  P = [x;y;z] <- update to approximate to solution
%           X(1) = [B(1) - {A(1,2)*P(2)+A(1,3)*P(3)+...+A(1,n)*P(n)}]/A(1,1)
%           X(2) = [B(2) - {A(2,1)*P(1)+A(1,3)*P(3)+...+A(1,n)*P(n)}]/A(2,2)
%                         ...
%           X(k) = [B(k) - {A(k,1)*P(1)+...+A(k,k-1)*P(k-1)
%                          +A(k,k+1)*P(k+1)+...+A(1,n)*P(n)}]/A(k,k)
%                         ...
%           X(n) = [B(n) - {A(n,1)*P(1)+...A(n,n-1)*P(n-1)}]/A(n,n)           
%
%--------------------------------------------------------------------------
%      matlab : 
%      >>   A=[4,-1,1; 4,-8,1;-2,1,5]
%      >>   B=[7;-21;15]
%      >>   P=[1;1;1]
%      >>   del = 0.00001
%      >>   max = 100
%      >>   X=jacobi(A,B,P,del,max)
%           => sol=(2,4,3)
%
%      In the below, vector form
%--------------------------------------------------------------------------
  function X=jacobi(A,B,P,del,max)
   
     N=length(B);
     for k=1:max
         for j=1:N
             X(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1;j+1:N]))/A(j,j);
         end
         err=norm(X'-P);
         relerr=err/((norm(X)+del))
         P=X';
         if (err<del)|(relerr<del)
            break
         end
     end
     X=X'
 
%--------------------------------------------------------------------
