qq_bicgstab.m

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

function [x,iter,flg,operation,res]=QQ_bicgstab(A,b,tol,maxit,u,q,x0)
% BI-CGSTAB algorithm
% [x,iter,flg,operation,res]=Q_bicgstab(A,b,tol,maxit,u,q,x0)
%  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        : 2006-07-08
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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(u)
    existu=1;
    nz_u=nnz(u);
    if ~isequal(size(u),[m,n])
        error('The preconditioners r should match the size of A!')
    end
else
    existu=0;
end
if nargin>=6 & ~isempty(q)
    existq=1;
    nz_q=nnz(q);
%     if ~isequal(size(q),[m,n])
%         error('The preconditioners D should match the size of A!')
%     end
else
    existq=0;
end

if nargin==7 & ~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>7
    error('Too many input arguments!')
end


% the main process of BICGSTAB
% the set up for this method
flg=1;
x=x0;
r=b-A*x;
r_w=r;
operation(1)=nz;iter=0;res(1)=1;
for i= 1: maxit
    operation(i+1)=operation(i);
    rho=r_w'*r;
    if rho==0
        flg=2; % method fails
        return;
    end
    if i==1
        p=r;
    else
        my_beta=(rho/rho1)*(my_alpha/w);
        p=r+my_beta*(p-w*v);
    end
    
    %preconditioning
    p_w=p;
    if existq
        p_w=q*p_w;
        operation(i+1)=operation(i+1)+nz_q;
    end
    if existu
        p_w=u\p_w;
        operation(i+1)=operation(i+1)+nz_u;
    end    
    v=A*p_w;
    my_alpha=rho/(r_w'*v);
    s=r-my_alpha*v;
    
    norm_s=norm(s);
    res(i+1)=norm_s/norm_b;
    iter=iter+0.5;
    if norm_s<tol % converges
        operation(i+1)=operation(i+1)+5*n+nz;
        flg=0;
        x=x+my_alpha*p_w;
        break;
    end
    % Precondition
    s_w=s;    
    if existq
        s_w=q*s_w;
        operation(i+1)=operation(i+1)+nz_q;
    end
    if existu
        s_w=u\s_w;
        operation(i+1)=operation(i+1)+nz_u;
    end
    t=A*s_w;
    w=t'*s/(t'*t);
    x=x+my_alpha*p_w+w*s_w;
    r=s-w*t;
    norm_r=norm(r,2);
    res(i+1)=norm_r/norm_b;
    operation(i+1)=operation(i+1)+10*n+2*nz;
    rho1=rho;
    iter=iter+0.5;
    if norm_r<tol  %converges
        flg=0;
        break;
    end
end
if flg~=0
    iter=maxit;
end