Shifted BiCG method

shifted_bicg.m

function x = shifted_bicg(A, b, sigma, tol, max_iter)

%%% initialization
if nargin < 4
    tol = 1e-6;
    max_iter = 1000;
end

if nargin < 5
    max_iter = 1000;
end

shift_convergence_flag = zeros(1, length(sigma));
iters = ones(1, length(sigma)) * max_iter;
norm_b = norm(b);

r = b;
r_shadow = r;

x_seed = zeros(length(A), 1);
x_shift = zeros(length(A), length(sigma));

p = zeros(length(A), 1);
p_shadow = zeros(length(A), 1);
p_shift = zeros(length(A), length(sigma));

alpha_old = 1;
rho_old = 1;
pi_shift_old = ones(1, length(sigma));
pi_shift = ones(1, length(sigma));

for iter = 1:max_iter
    %%% seed system
    rho = r_shadow' * r;
    beta = -rho / rho_old;
    p = r - beta * p;
    p_shadow = r_shadow - conj(beta) * p_shadow;
    q = A * p;
    alpha = rho / (r_shadow' * q);
    x_seed = x_seed + alpha * p;
    
    %%% shifted system
    pi_shift_new = (1 + alpha * sigma) .* pi_shift ...
                + alpha * beta / alpha_old * (pi_shift_old - pi_shift);
    beta_shift = power(pi_shift_old ./ pi_shift, 2) * beta;
    alpha_shift = pi_shift ./ pi_shift_new * alpha;
    for i = find(shift_convergence_flag == 0)
        p_shift(:,i) = r / pi_shift(i) - p_shift(:,i) * beta_shift(i);
        x_shift(:,i) = x_shift(:,i) + p_shift(:,i) * alpha_shift(i);
    end
    
    %%% update residual and scalars
    r = r - alpha * q;
    r_shadow = r_shadow - conj(alpha) * A' * p_shadow;
    alpha_old = alpha;
    rho_old = rho;
    pi_shift_old = pi_shift;
    pi_shift = pi_shift_new;
    
    for i = find(shift_convergence_flag == 0)
        if norm(r/pi_shift(i))/norm_b < tol
            shift_convergence_flag(i) = 1;
            iters(i) = iter;
        end
    end
    
    %%% とりあえずseedが収束したら終わりにする
    if norm(r) / norm_b < tol
        break
    end
    
end

x = [x_seed x_shift];

%%% display results
disp('results for seed systems...');
disp(strcat('relative residual:', num2str(norm(r)/norm(b))));
disp(strcat('true residual:', num2str(norm(b-A*x_seed)/norm(b))));

disp('results for shifted systems...');
for i = 1:length(sigma)
    disp(strcat('sigma:', num2str(sigma(i))));
    disp(strcat('relative residual:', num2str(norm(r/pi_shift(i))/norm(b))));
    disp(strcat('true residual:', num2str(norm(b-(A+sigma(i)*eye(length(A)))*x_shift(:,i))/norm(b))));
end