function x=pcg(errtol)%(A,b,x,omiga)
%
omiga=1;
tic;   %计时器开始，与计算无关
[A,b]=prosparse; %生成矩阵A,b 即函数为Ax=b

x=ones(100,1);% x的初始值
[D,L,U]=dlu(A);% 对A进行分解子函数，见下面
M=(D-omiga*L)*inv(D)*(D-omiga*U)/(omiga*(2-omiga));%计算ssor预处理矩阵M

%下面为pcg方法的计算过程，完全翻译，所以不做标注
r=b-A*x;
invm=inv(M);
z=invm*r;
p=z;
x0=x-1;
k=1;
while (k<=50)&&(max(abs(x-x0))>errtol)
    zr=z'*r;
    alpha=zr/(p'*(A*p));
    x0=x;
    x=x+alpha*p;
    r=r-alpha*A*p;
    z=invm*r;
    beta=(z'*r)/zr;
    p=z+beta*p;
    k=k+1;
end
disp(['k=',num2str(k)]);%先是迭代次数k和耗时
disp(['pcg method (s) :',num2str(toc)]);

function [D,L,U]=dlu(A) %dlu分解子函数
D=zeros(size(A));
L=D;
U=D;

for i=1:size(A,1)
    D(i,i)=A(i,i);
end
for i=1:size(A,1)
    for j=1:size(A,2)
        L(i,j)=-(i<j)*A(i,j);
        U(i,j)=-(i>j)*A(i,j);
    end
end

function [A,f]=prosparse  %生成待机算数组子函数，即模型的离散化

Aii=eye(10)*4-diag(ones(1,9),1)-diag(ones(1,9),-1);
A=zeros(100);
A(1:10,1:20)=[Aii,-eye(10)];
for i=11:10:81
    A(i:i+9,i-10:i+19)=[-eye(10),Aii,-eye(10)];
end
A(91:end,81:end)=[-eye(10),Aii];
A=sparse(A);  %把矩阵A转化为系数矩阵以减小存储量，可以不用

f=zeros(10);%生成函数右端部分f。
for i=1:10
    for j=1:10
        f(i,j)=(i^2+j^2)*exp(i*j/121)/(11^4);
    end
end
f(1,:)=f(1,:)+1;%这四行用来进行边界处理
f(:,1)=f(:,1)+1;
f(end,:)=f(end,:)+exp((1:10)/11);
f(:,end)=f(:,end)+exp((1:10)'/11);
f=f(:);