my_bicg.m

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

function [x,iter,flg]=my_BiCG(A,b,tol,maxit,M1,D,M2,x0)
% Algorithm: BICG from Templates
% 
% [x,iter,flg]=my_BiCG(A,b,tol,maxit,M1,D,M2,x0)
% The input arguments are the usual ones
% flg=0 success; flg=1 fail without convergence; 
% flg=2 fail due to breakdown
% 
% Coding: Plum_LiLiang (SAM of UESTC)
% Date  : 2006-05-28 10:00 pm.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

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;
    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;
    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;
    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>8
    error('Too many input arguments!')
end

% The main process of iteration (BiCG)
% The set up
x=x0;
r=b-A*x0;
flg=1;
iter=maxit;
r_T=r;
for i=1:maxit
     % Preconditioning
    if existM1
        y=M1\r;
        y_T=M1'\r_T;
    else
        y=r;
        y_T=r_T;
    end
    if existD
        y=D\y;
        y_T=D'\y_T;
    end
    if existM2
        z=M2\y;
        z_T=M2'\y_T;
    else
        z=y;
        z_T=y_T;
    end
    rho1=z'*r_T; % The concept of BiConjugate
    if rho1==0   % break down
        flg=2;
        break;
    end
    if i==1
        p=z;
        p_T=z_T;
    else
        beta=rho1/rho;
        p=z+beta*p;
        p_T=z_T+beta*p_T;
    end
    q=A*p;
    q_T=A'*p_T;
    alpha=rho1/(p_T'*q); % The concept of BiConjugate
    x=x+alpha*p;
    r=r-alpha*q;
    r_T=r_T-alpha*q_T;
    norm_r=norm(r);
    if norm_r<=tol % check convergence
        flg=0;
        iter=i;
        break;
    end
    rho=rho1;
end