% This is matlab code of algorithm 9.2. Last updated on 6-17-2011.

function [x,err,iter,flag] = mlbicgstab(A,x,b,Q,M,max_it,tol,kappa)

% input: A: N-by-N matrix
%        M: N-by-N preconditioner matrix
%        Q: N-by-n auxiliary matrix with columns q_1,..., q_n
%        x: initial guess
%        b: right hand side vector
%        max_it: maximum number of iterations
%        tol: error tolerance
%        kappa: (real number) minimization step controller:
%               zero, standard minimization
%               positive, Sleijpen-van der Vorst minimization
% output: x: solution computed
%         err: error norm
%         iter: number of iterations performed
%         flag: = 0, solution found to tolerance
%               = 1, no convergence given max_it
%               = -1, breakdown
% storage: c: 1-by-n matrix
%          x, r, b, u_t, z: N-by-1 matrices
%          A, M: N-by-N matrices
%          Q, G, W: N-by-n matrices

N = size(A,2);
n = size(Q,2);
G = zeros(N,n);  % initialize work spaces
W = zeros(N,n);  
c = zeros(1,n);  % end initialization

iter = 0;
flag = 1;
bnrm2 = norm(b);
if bnrm2 == 0.0
    bnrm2 = 1.0;
end
r = b - A*x;
err = norm( r ) / bnrm2;
if err < tol
    flag = 0;
    return
end

G(:,1) = M \ r;
W(:,1) = A*G(:,1);
c(1) = Q(:,1)'*W(:,1);
if c(1) == 0
    flag = -1;
    return
end
e = Q(:,1)'*r;

for j = 0:max_it
    for i = 1:n-1
        alpha = e / c(i);
        x = x + alpha*G(:, i);
        r = r - alpha*W(:, i);
        err = norm(r)/bnrm2;
        iter = iter + 1;
        if err < tol
            flag = 0;
            return
        end
        if iter >= max_it
            return
        end

        e = Q(:,i+1)'*r;
        if j >= 1
            beta = -e / c(i+1);
            W(:,i+1) = r + beta*W(:,i+1);
            G(:,i+1) = beta*G(:,i+1);
            for s = i+1 : n-1
                beta = -Q(:,s+1)'*W(:,i+1)/c(s+1);
                W(:,i+1) = W(:,i+1) + beta*W(:,s+1);
                G(:,i+1) = G(:,i+1) + beta*G(:,s+1);
            end
            G(:,i+1) = (M \ W(:,i+1)) - (1/omega)*G(:,i+1);
        else
            G(:,i+1) = M \ r;
        end
        W(:,i+1) = A*G(:,i+1);
        for s = 0:i-1
            beta = -Q(:,s+1)'*W(:,i+1) / c(s+1);
            W(:,i+1) = W(:,i+1) + beta*W(:,s+1);
            G(:,i+1) = G(:,i+1) + beta*G(:,s+1);
        end
        c(i+1) = Q(:,i+1)'*W(:,i+1);
        if c(i+1) == 0
            flag = -1;
            return
        end
    end
    alpha = e / c(n);
    x = x + alpha*G(:, n);
    r = r - alpha*W(:,n);
    err = norm(r)/bnrm2;
    if err < tol
        flag = 0;
        iter = iter + 1;
        return
    end
    u_t = M \ r;
    z = A*u_t;
    omega = z'*z;
    if omega == 0
        flag = -1;
        return
    end
    rho = z'*r;
    omega = rho / omega;
    if kappa > 0
        rho = rho / (norm(z)*norm(r));
        abs_rho = abs(rho);
        if (abs_rho < kappa) & (abs_rho ~= 0)
            omega = omega*kappa/abs_rho;
        end
    end
    if omega == 0
        flag = -1;
        return
    end
    x = x + omega*u_t;
    r = r - omega*z;
    err = norm(r)/bnrm2;
    iter = iter + 1;
    if err < tol
        flag = 0;
        return
    end
    if iter >= max_it
        return
    end
    e = Q(:,1)'*r;
    beta = - e/c(1);
    W(:,1) = r + beta*W(:,1);
    G(:,1) = beta*G(:,1);
    for s = 1:n-1
        beta = -Q(:,s+1)'*W(:,1)/c(s+1);
        W(:,1) = W(:,1) + beta*W(:,s+1);
        G(:,1) = G(:,1) + beta*G(:,s+1);
    end
    G(:,1) = (M \ W(:,1)) - (1/omega)*G(:,1);
    W(:,1) = A*G(:,1);
    c(1) = Q(:,1)'*W(:,1);
    if c(1) == 0
        flag = -1;
        return
    end
end