clear all;
format long;  %修改默認數據類型
tic;   %開始計時
n=1000;   %矩陣規模
l=500; %預處理子矩陣的個數
e=0.00000001;  %精度
A=hilb(n);   %生成系數矩陣
K=cond(A);  %條件數
M=zeros(n);
for i=1:1:l
   for j=1:1:(n/l)
       for k=1:1:(n/l)
           line=j+(i-1)*(n/l);
           low=k+(i-1)*(n/l);
           M(line,low)=A(line,low);
       end
   end
end
x=ones(n,1);   %生成精確解
b=A*x;   %生成精確解對應的b
X=zeros(n,1);  %初始
ek1=sqrt((X-x)'*A*(X-x)); %開始是不會知道精確解的！
sulv=0;  % 評價指標
r=b-A*X;
ni=inv(M);
z=ni*r;
p=z;
k=1;%迭代次數
fan=norm(r);
while (k<=1000)&&(fan>e)
zr=r'*z;
afa=zr/(p'*(A*p));
X=X+afa*p;
ek2=sqrt((X-x)'*A*(X-x));
sulv=sulv+ek2/ek1;
ek1=ek2;
r=r-afa*(A*p);
z=ni*r;
fan=norm(r);
beta=r'*z/zr;
p=z+beta*p;
k=k+1;
end
en=norm(X-x)  %求解的誤差
sulv=sulv/(k-1)  %下降速率的平均值
B=A*X;
toc;