my_bicgstab.m

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

function [x,iter,flg,operation,res]=my_bicgstab(A,b,tol,maxit,M1,D,M2,x0,p_type)
% BI-CGSTAB algorithm
% [x,iter,flg,operation,res]=my_bicgstab(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        : 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(M1)
    existM1=1;
    nz_M1=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;
    nz_D=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;
    nz_M2=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 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;i=0;err=1;
while err>tol
    i=i+1;
    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 p_type==0  % implicit preconditioner
        if existM1
            p_w=M1\p_w;
            operation(i+1)=operation(i+1)+nz_M1;
        end
        if existD
            p_w=D\p_w;
            operation(i+1)=operation(i+1)+nz_D;
        end
        if existM2
            p_w=M2\p_w;
            operation(i+1)=operation(i+1)+nz_M2;
        end
    elseif p_type==1 % explicit preconditioner
        if existM1
            p_w=M1*p_w;
            operation(i+1)=operation(i+1)+nz_M1;
        end
        if existD
            p_w=D*p_w;
            operation(i+1)=operation(i+1)+nz_D;
        end
        if existM2
            p_w=M2*p_w;
            operation(i+1)=operation(i+1)+nz_M2;
        end
    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 p_type==0  % implicit preconditioner
        if existM1
            s_w=M1\s_w;
            operation(i+1)=operation(i+1)+nz_M1;
        end
        if existD
            s_w=D\s_w;
            operation(i+1)=operation(i+1)+nz_D;
        end
        if existM2
            s_w=M2\s_w;
            operation(i+1)=operation(i+1)+nz_M2;
        end
    elseif p_type==1 % explicit preconditioner
        if existM1
            s_w=M1*s_w;
            operation(i+1)=operation(i+1)+nz_M1;
        end
        if existD
            s_w=D*s_w;
            operation(i+1)=operation(i+1)+nz_D;
        end
        if existM2
            s_w=M2*s_w;
            operation(i+1)=operation(i+1)+nz_M2;
        end
    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
    err=norm_r;
    if iter==maxit
        flg=1;
        break;
    end
end
res = res/norm_b;
% if flg~=0
%     iter=maxit;
% end