cg.m

求解线性系统的Krylov方法的工具箱

%function y=CG(A,b,x0)
%developped by Plum_Liliang
%this algorithm is to sovle the linear system Ax=b
%where A is a n-by-n matrix
function [y,x]=CG(A,b,x0)
k=0;n=length(A);
if nargin<3
    x0=zeros(n,1);
end
r0=b-A*x0;
r1=r0;
err=1;
while err>0.00001
    k=k+1;
    if k<100
      if k==1
          p=r0;
      else
%           beta=-p'*A*r0/(p'*A*p);
%           p=r0+beta*p;
          beta=r1'*r1/(r0'*r0);
          p=r0+beta*p;
      end  
%       alpha=p'*r0/(p'*A*p);
%       x1=x0+alpha*p;
                alpha=p'*r1/(p'*A*p);
                x1=x0+alpha*p;
      r0=r1;
      r1=r1-alpha*A*p;
      err=norm(r1,inf);
      x0=x1;
    else
        break;
    end
end
y=k;
x=x1;