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

% mlbicgstab.m solves the linear system Ax=b using the
% multiple Lanczos BiCGSTAB Method with preconditioning.
%
% input   A        matrix
%         M        preconditioner matrix
%         Q        initial matrix [q_1,...,q_k] with columns q_1,...,q_k
%         x        initial guess vector
%         b        right hand side vector
%         max_it   (integer) maximum number of iterations
%         tol      error tolerance
%
% output  x        solution vector
%         error    error norm
%         iter     (integer) number of iterations performed
%         flag     (integer): 0 = solution found to tolerance
%                             1 = no convergence given max_it
%                            -1 = breakdown

  n = size(A,2);
  k = size(Q,2);
  iter = 0;
  flag = 1;
  jmax_it = fix(max_it./k);
  G = zeros(n,k);
  W = zeros(n,k);
  D = zeros(n,k);

  bnrm2 = norm(b);
  if bnrm2 == 0.0
    bnrm2 = 1.0; 
  end
  r = b - A*x;
  G(:,k) = r;

  error = norm( r ) / bnrm2;
  if error < tol 
    flag = 0; 
    return
  end

  for j = 0:jmax_it

    if (iter >= max_it) | (flag == -1) | (flag == 0)
      break
    end

    g_tld = M\G(:,k);
    W(:,k) = A*g_tld;
    c(k) = Q(:,1)'*W(:,k);

    if c(k) == 0 
      flag = -1;
      break
    end
    alpha = Q(:,1)'*r/c(k);

    u = r - alpha*W(:,k);
    u_tld = M\u;
    temp = A*u_tld;

    rho =temp'*temp;
    if rho == 0
      flag = -1;
      break
    end
    rho = -u'*temp/rho;

    x = x - rho*u_tld + alpha*g_tld;
    iter = iter + 1;
    r = rho*temp + u;
    error = norm(r)/bnrm2;
    
    if error < tol
      flag = 0;
      break
    end

    for i = 1:k
      
      if iter >= max_it
        break
      end

      zd = u;
      zg = r;
      zw = zeros(n,1);

      if j ~= 0
        for s = i:k-1
          beta = -Q(:,s+1)'*zd/c(s);
          zd = zd + beta*D(:,s);
          zg = zg + beta*G(:,s);
          zw = zw + beta*W(:,s);
        end
      end

      beta = rho*c(k);
      if beta == 0
        flag = -1;
        break
      end
      beta = -Q(:,1)'*(r + rho*zw)/beta;
      zg = zg + beta*G(:,k);
      zw = rho*(zw + beta*W(:,k));
      zd = r + zw;

      for s = 1:i-1
        beta = -Q(:,s+1)'*zd/c(s);
        zd = zd + beta*D(:,s);
        zg = zg + beta*G(:,s);
      end

      D(:,i) = zd - u;
      G(:,i) = zg + zw;

      if i < k
        c(i) = Q(:,i+1)'*D(:,i);
        if c(i) == 0
          flag = -1;
          break
        end
        alpha = Q(:,i+1)'*u/c(i);
        u = u - alpha*D(:,i);
        g_tld = M\G(:,i);

        x = x + rho*alpha*g_tld;
        iter = iter + 1;
        W(:,i) = A*g_tld;
        r = r - rho*alpha*W(:,i);
        error = norm(r)/bnrm2;

        if error < tol
          flag = 0;
          break
        end
      end
    end
  end