cg.m

来自「％Jacobi迭代法求方程组等程序」· M 代码 · 共 18 行

function x=CG(A,x0,b)
%CG算法
%A=[6,3;3,2];b=[0,-1]';x0=[0,0]';tc_CG(A,x0,b)
[m,n]=size(A);
if m~=n
    fprintf('\n"A" 应该为一个矩阵！')
end   %if
while(1)
r0=b-A*x0;p0=r0;
a0=r0'*r0/((A*p0)'*p0);
x=x0+a0*p0;
r=r0-a0*A*p0;
  if  (r'*r)^0.5<=1e-10
    return  
  end   %if
p0=r+(r'*r/(r0'*r0))*p0;
x0=x;r0=r;
end     %while