谁能帮我编一个线性方程Ax=b的共轭梯度法的程序？谢谢了

function [x,n]= conjgrad (A,b,x0) %共轭梯度法求线性方程组Ax=b的解

r1 = b-A*x0

p = r1

n = 0

for i=1:rank(A) %以下过程可参考算法流程

if(dot(p,A*p) <1.0e-50) %循环结束条件

break

end

alpha = dot(r1,r1)/dot(p,A*p)

x = x0+ alpha*p

r2 = r1- alpha*A*p

if(r2 <1.0e-50)%循环结束条件

break

end

belta = dot(r2,r2)/dot(r1,r1)

p = r2+belta*p

n = n + 1

end

function [x,n]= preconjgrad (A,b,x0,M,eps) %预处理共轭梯度法求线性方程组Ax=b的解

if nargin == 4

eps = 1.0e-6

end

r1 = b-A*x0

iM = inv(M)

z1 = iM*r1

p = z1

n = 0

tol= 1

while tol>=eps

alpha = dot(r1,z1)/dot(p,A*p)

x = x0 + alpha*p

r2 = r1 - alpha*A*p

z2 = iM*r2

belta = dot(r2,z2)/dot(r1,z1)

p = z2+belta*p

n = n + 1

tol = norm(x-x0)

x0 = x %更新迭代值

r1 = r2

z1 = z2

end

function A=xishu(N) %存储离散化后非边界点构成的系数矩阵

A=zeros(N^2)

for i=1:N^2

A(i,i)=4

if i+N<N^2+1

A(i,i+N)=-1

A(i+N,i)=A(i,i+N)

end

if mod(i,N)~=0

A(i,i+1)=-1

A(i+1,i)=A(i,i+1)

end

end

function x=cg(a,b,x)

%共轭向量法求解线性方程组

% 初始输入变量 a：系数矩阵 b：方程右端向量 x 输入的初始解

n=length(b)

r=b-a*x

p=r

q0=r'*r

while q0>1e-10

w=a*p

t=q0/(p'*w)

x=x+t*p

r=r-t*w

q=r'*r

s=q/q0

p=r+s*p

q0=q

end

---