my_cgnr.m

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

function [x,iter,flg,res,operation]=my_CGNR(A,b,tol,maxit,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);
flg=0;
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;rw=A'*r;
res(1)=1;
operation(1)=nz+nz;
iter=1;
% preconditioning
z=rw;
if ptype==0 % implicit precondition
    if existM1
        z = M1\z;
        operation(1)=operation(1)+nzM1;
    end
    if existD
        z = D\z;
        operation(1)=operation(1)+nzD;
    end
    if existM2
        z = M2\z;
        operation(1)=operation(1)+nzM2;
    end
%     explicit precondition
elseif ptype==1
    if existM1
        z = M1*z;
        operation(1)=operation(1)+nzM1;
    end
    if existD
        z = D*z;
        operation(1)=operation(1)+nzD;
    end
    if existM2
        z = M2*z;
        operation(1)=operation(1)+nzM2;
    end
end
p=z;err=1;zr1=z'*rw;
operation(1)=operation(1)+n;
while err>tol
    iter=iter+1;
    w=A*p;norm_w=w'*w;
    alp=zr1/norm_w;
    x=x+alp*p;
    r=r-alp*w;
    err=norm(r)/norm_b;
    res(iter)=err;
    rw=A'*r;
    operation(iter)=2*nz+operation(iter-1)+2*n;
    z=rw;
    if ptype==0 % implicit precondition
        if existM1
            z = M1\z;
            operation(iter)=operation(iter)+nzM1;
        end
        if existD
            z = D\z;
            operation(iter)=operation(iter)+nzD;
        end
        if existM2
            z = M2\z;
            operation(iter)=operation(iter)+nzM2;
        end
        %     explicit precondition
    elseif ptype==1
        if existM1
            z = M1*z;
            operation(iter)=operation(iter)+nzM1;
        end
        if existD
            z = D*z;
            operation(iter)=operation(iter)+nzD;
        end
        if existM2
            z = M2*z;
            operation(iter)=operation(iter)+nzM2;
        end
    end
    zr2=z'*rw;beta=zr2/zr1;zr1=zr2;
    p=z+beta*p;
    if iter>maxit
        flg=1;
        break;
    end
end
iter=iter-1;