clear all;
format long;  %修改默認數據類型
tic;   %開始計時
n=1000;   %矩陣規模
e=0.00000001;  %精度
A=hilb(n);   %生成系數矩陣
K=cond(A);  %條件數
M=zeros(n);
sum=0;
for i=1:1:45  %一共45個子矩陣
   for j=1:1:i
       for k=1:1:i
           line=j+sum;
           low=k+sum;            
           if(line<=n)&&(low<=n)
               M(line,low)=A(line,low);
           end
       end
   end
   sum=sum+i;
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<=100)&&(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;