my_rpsqmr.m

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

function [x,iter,flg,operation,res]=my_RPSQMR(A,b,tol,maxit,M1,D,M2,x0,p_type)
% simplyfied QMR algorithm
% [x,iter,flg,operation,res]=my_SQMR(A,b,tol,maxit,M1,D,M2,x0,p_type)
%  flg=0 success; flg=1 not convergent; flg=2 method fails
% p_type=0 implicit preconditioner
% p_type=1 explicit preconditioner

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

% check the input arguments
if nargin<2
    error('Not enough input arguments!');
end

[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

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

norm_b=norm(b,2);
tol=tol*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

% Check the initial guess x0
if nargin==8 & ~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 | isempty(p_type)
    p_type=0;
end

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


% the main process of SQMR
% the set up for this method
x=x0;
r=b-A*x;norm_r0=norm(r);
operation(1)=nz;
iter=1;
q=r;
flg=1;d=zeros(n,1);
% Precondition
if p_type==0 % implicit precondition
    if existM1
        q = M1\q;
        operation(1)=operation(1)+nzM1;
    end
    if existD
        q = D*q;
        operation(1)=operation(1)+nzD;
    end
    if existM2
        q = M2\q;
        operation(1)=operation(1)+nzM2;
    end
%     explicit precondition
elseif p_type==1
    if existM1
        q = M1*q;
        operation(1)=operation(1)+nzM1;
    end
    if existD
        q = D*q;
        operation(1)=operation(1)+nzD;
    end
    if existM2
        q = M2*q;
        operation(1)=operation(1)+nzM2;
    end
end
tau=norm(r);
rho=r'*q;
if isinf(rho) | isnan(rho)
    iter=0;
    flg=2;
    res=0;
    operation=0;
    return;
end
err=norm_b;res(1)=norm_b;
v1=0;
while err>tol
    iter=iter+1;
    t=A*q; delta=q'*t;
%     res(iter)=delta;
    operation(iter)=operation(iter-1)+nz+n;
%     err=delta;
%     if abs(delta)<tol
%         res(iter)=abs(delta);
%         flg=0;
%         iter=iter-1;
%         res=res/res(1);
%         return;
%     end
    alf=rho/delta;
    r=r-alf*t;
    norm_r=norm(r,2);
    err=norm_r;
    res(iter)=norm_r;
    if abs(norm_r)<tol
        flg=0;
        iter=iter-1;
        res=res/res(1);
        return;
    end
    v=norm(r)/tau;
    c=1/sqrt(1+v);
    tau=tau*v*c;
    d=(c^2)*(v1^2)*d+(c^2)*alf*q;
    x=x+d;
    v1=v;
%     if abs(rho)<tol
%         res(iter)=abs(rho);
%         flg=0;
%         iter=iter-1;
%         res=res/res(1);
%         return;
%     end
    u=r;
    % Precondition
    if p_type==0 % implicit precondition
        if existM1
            u = M1\u;
            operation(iter)=operation(iter)+nzM1;
        end
        if existD
            u = D*u;
            operation(iter)=operation(iter)+nzD;
        end
        if existM2
            u = M2\u;
            operation(iter)=operation(iter)+nzM2;
        end
        %     explicit precondition
    elseif p_type==1
        if existM1
            u = M1*u;
            operation(iter)=operation(iter)+nzM1;
        end
        if existD
            u = D*u;
            operation(iter)=operation(iter)+nzD;
        end
        if existM2
            u = M2*u;
            operation(iter)=operation(iter)+nzM2;
        end
    end
    rho1=r'*u;beta=rho1/rho;rho=rho1;
    operation(iter)=operation(iter)+2*n;
    q=u+beta*q;
    if iter>=maxit
        flg=1;
        break;
    end
end