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)

break;

end

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

x = x0+ alpha*p;

r2 = r1- alpha*A*p;

if(r2

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