my_csym.m

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

function [x,iter,flg,res]=my_csym(A,b,tol,maxit,M1,D,M2,x0)
% This code implements CSYM algorithm LAA 287(1999) 105-123
% with preconditioners M1 and M2, D if necessary
% flg=0 success; flg=1 failed to convergent whin maxit iterations;
% flg = 2 breakdown.

% Developed by: Liang Li, UESTC China
% Date        : 2008-01-30

% 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
nz=nnz(A);
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;
    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

% set up of CSYM
x = x0;
flg = 1;
r = b - A*x;
norm_r0 = norm(r);
tolabs = tol*norm_r0;
q0 = zeros(n,1); q1 = conj(r)/norm_r0;
tempv = A*q1;
if existM1
    tempv = M1 \ tempv;
end
if existD
    tempv = D \ tempv;
end  
if existM2
    tempv = M2 \ tempv; 
end
alf1 = q1.'*tempv; beta1 = 0;
c_1 = 0; s_1 = 0; c0 = 1; s0 = 0;
p0 = zeros(n,1);
p_1 = p0;
tau1 = norm_r0;
iter = maxit;
% The main loop of csym
for k = 1:maxit
    eta = c_1*c0*beta1 + conj(s0)*alf1;
    theta = conj(s_1)*beta1;
    gamma = c0*alf1 - c_1*s0*beta1;
    w = tempv - alf1*conj(q1) - beta1*conj(q0);
    beta1 = norm(w);
    if beta1 == 0
        flg = 2;
        iter = k;
        return;
    end
    q0 = q1;
    q1 = conj(w)/beta1;
    tempv = A*q1;
    if existM1
        tempv = M1 \ tempv;
    end
    if existD
        tempv = D \ tempv;
    end  
    if existM2
        tempv = M2 \ tempv; 
    end
    if gamma ~= 0
        deta = sqrt(abs(gamma)^2+beta1^2);
        c = abs(gamma)/deta;
        s = conj(gamma)*beta1/(abs(gamma)*deta);
        xi = gamma*deta/abs(gamma);
    else
        c = 0; s = 1; xi = beta1;
    end    
    p1 = (q0 - eta*p0 - theta*p_1)/xi;
    p_1 = p0; p0 = p1;
    x = x + tau1*c*p1;
    tau1 = -s*tau1;
    res(k) = abs(tau1);
    if res(k)<tolabs
        flg = 0;
        iter = k;
        return;
    end
end