function [x,N]= BJ(A,b,x0,d,eps,M)
if nargin==4
    eps= 1.0e-6;
    M  = 10000;
elseif nargin<4
    error
    return
elseif nargin ==5
    M  = 10000;             %參數的默認值
end  

NS = size(A);
n = NS(1,1);                
if(sum(d) ~= n)
    disp('分塊錯誤！');
    return;
end
bnum = length(d);
bs = ones(bnum,1);
for i=1:(bnum-1)
    bs(i+1,1)=sum(d(1:i))+1;  
    %獲得對角線上每個分塊矩陣元素索引的起始值
end

DB = zeros(n,n);                
for i=1:bnum
    endb = bs(i,1)+d(i,1)-1;
    DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb);
    %求A的對角分塊矩陣
end
for i=1:bnum
    endb = bs(i,1)+d(i,1)-1;
    invDB(bs(i,1):endb,bs(i,1):endb)=inv(DB(bs(i,1):endb,bs(i,1):endb));
    %求A的對角分塊矩陣的逆矩陣
end

N = 0;
tol = 1;
while tol>=eps 
    x = invDB*(DB-A)*x0+invDB*b;    %由于LB+DB=DB-A
    N = N+1;                               %迭代步數
    tol = norm(x-x0);                   %前后兩步迭代結果的誤差
    x0 = x;
    if(N>=M)
        disp('Warning: 迭代次數太多，可能不收斂！');
        return;
    end
end