my_cgnr.asv

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

function [x,iter,flg,res,operation]=my_CGNR(A,b,tol,M1,D,M2,ptype,x0)
% [x,iter,flg,res,operation]=my_CGNR(A,b,tol,M1,D,Mr2,ptype,x0)
% This code implements CG algorithm with preconditioners M1 and M2, D if
% necessary
% flg=0 success; flg=1 fail;
% res--a vector record the residuals for each step
% operation--a vector record the number of operations for each step

% Developed by: Plum_Liliang UESTC China
% Date        : 2007-01-12
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check the dimension of A and the right-hand vector b
[m,n]=size(A);
if m~=n
    error('The coefficient matrix must be square!')
end
[m_b,n_b]=size(b);
if n_b~=1
    error('b must be a vector!')
end
nz=nnz(A);
if m~=m_b
    error('The right-hand vecotor must have the same length of A!')
end

% Check the input arguments and assign the default value
if nargin<2
    error('Not enough input arguments!');
end

if (nargin < 3) | isempty(tol)
   tol = 1e-6;
end
if (nargin < 4) | isempty(maxit)
   maxit = min(n,20);
end

% if 'b' is zero then the solution is zero
norm_b=norm(b);
if norm_b==0
    x=0;
    iter=0;
    x=0;
    return;
end

% Check the preconditioners
if nargin>=5 & ~isempty(M1)
    existM1=1;nzM1=nnz(M1);
    if ~isequal(size(M1),[m,n])
        error('The preconditioners M1 should match the size of A!')
    end
else
    existM1=0;
end
if nargin>=6 & ~isempty(D)
    existD=1;nzd=nnz(D);
    if ~isequal(size(D),[m,n])
        error('The preconditioners D should match the size of A!')
    end
else
    existD=0;
end
if nargin>=7 & ~isempty(M2)
    existM2=1;nzM2=nnz(M2);
    if ~isequal(size(M2),[m,n])
        error('The preconditioners M2 should match the size of A!')
    end
else
    existM2=0;
end

if nargin<8 &~isempty(ptype)
    ptype=0; % default preconditioning type is the usual case
end

% Check the initial guess x0
if nargin==9 & ~isempty(x0)
    if ~isequal(size(x0),[m,1])
        error('The initial guess x0 must have the same length of A')
    end
else
    x0=zeros(m,1);
end

if nargin>9
    error('Too many input arguments!')
end

% set up for CGNR
x=x0;
r=b-A*x;r=A'*r;
res(1)
% preconditioning

while err>tol
    k=k+1;
    if k<100
      if k==1
          p=A'*r;
      else
          ATr=A'*r;
          beta=((A'*r)'*(A'*r))/((A'*r0)'*(A'*r0));
          p=ATr+beta*p;
      end
      Ap=A*p;
      alpha=((ATr)'*(ATr))/((Ap)'*(Ap));
      x1=x0+alpha*p;
      r0=r;
      r=r-alpha*Ap;
      err=norm(r,inf);
      x0=x1;
    else
        break;
    end
end