my_chebyshev.m

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

function [x,iter,flg,operations]=my_chebyshev(A,b,tol,maxit,M1,D,M2,x0,p_type,lammax,lammin)
% Chebyshev semi-iterative algorithm
% [x,iter,flg,operations]=my_chebyshev(a,b,tol,maxit,M1,M2,x0,p_type,lammax,lammin)
% lammax: the maximum eigenvalue; lammin: the minimum one
% p_type=0 implicit preconditioner
% p_type=1 explicit preconditioner
%  flg=0 success; flg=1 not convergent

% Developed by: Plum_Liliang UESTC China
% Date        : 2006-07-09
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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);
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<10 | isempty(lammax)
   lammax=2; 
end

if nargin<11 | isempty(lammin)
    lammin=1e-3;
end

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


% setup of this method
flg=1;
d=(lammax+lammin)/2;
c=(lammax-lammin)/2;
r=b-A*x0;
operations=nz;
iter=0;
x=x0;
for i=1:maxit
    % Precondition
    if p_type==0
        if existM1==1
            r=M1\r;
            operations=operations+nz_M1;
        end
        if existD==1
            r=D\r;
            operations=operations+nz_D;
        end
        if existM2==1
            r=M2\r;
            operations=operations+nz_M2;
        end
    elseif p_type==1
        if existM1==1
            r=M1*r;
            operations=operations+nz_M1;
        end
        if existD==1
            r=D*r;
            operations=operations+nz_D;
        end
        if existM2==2
            r=M2*r;
            operations=operations+nz_M2;
        end
    end
    if i==1
        p=r;
        alpha1=2/d;
    else
        beta1=(c*alpha1/2)^2;
        alpha1=1/(d-beta1);
        p=r+beta1*p;
    end
    x=x+alpha1*p;
    r=b-A*x;
    norm_r=norm(r);
    operations=operations+2*n+nz;
    if norm_r<tol % convergent
        flg=0;
        iter=i;
        break;
    end
end
if flg==1 % not convergent
    iter=maxit;
end