function [x,err,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        N-by-N matrix
%         M        N-by-N preconditioner matrix
%         Q        N-by-n auxiliary matrix [q_1,...,q_n] with columns 
%                  q_1,...,q_n
%         x        initial guess 
%         b        right hand side vector
%         max_it   (integer) maximum number of iterations
%         tol      error tolerance
%
% output  x        solution 
%         err      error norm
%         iter     (integer) number of iterations performed
%         flag     (integer): 0 = solution found to tolerance
%                             1 = no convergence given max_it
%                            -1 = breakdown
%
% storage:  
%      alphaa, betaa, rhoo, bnrm2, err, tol, d, e, f ---  real/complex numbers
%      i, j, s, n, N, flag, iter, max_it ---  integers
%                                      c ---  real/complex array of length n
%      x, r, g, g_tld, u, u_tld, w, z, b ---  vectors of length N
%                                   A, M ---  N-by-N matrices
%                                Q, G, W ---  N-by-n matrices
%                                      D ---  N-by-(n-1) matrix. 
%                                             D is not required when n = 1.   
%

  N = size(A,2);
  n = size(Q,2);

  G = zeros(N,n);          % reserve spaces for G, W, D and c
  W = zeros(N,n);
  if n > 1
     D = zeros(N,n-1);
  end
  c = zeros(1,n);          % end reservation

  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(:,n) = r;
  g_tld = M\r;
  W(:,n) = A*g_tld;
  c(n) = Q(:,1)'*W(:,n);

  if c(n) == 0
      flag = -1;
      return
  end
  
  e = Q(:,1)'*r;

  for j = 0:max_it
    alphaa = e/c(n);
    u = r - alphaa*W(:,n);
    u_tld = M\u;
    z = A*u_tld;

    rhoo = z'*z;
    if rhoo == 0
      flag = -1;
      return
    end

    rhoo = - z'*u / rhoo;
    x = x - rhoo*u_tld + alphaa*g_tld;
    r = rhoo*z + u;

    err = norm(r)/bnrm2;
    iter = iter + 1;

    if err < tol
      flag = 0;
      return
    end

    if iter >= max_it
       return
    end

    d = rhoo*c(n);
    if d == 0
       flag = -1;
       return
    end

    for i = 1:n-1
      z = u;
      g = zeros(N,1);     
      w = zeros(N,1);
      f = Q(:,i+1)'*u;

      if j >= 1
         betaa = - f / c(i);
         z = z + betaa*D(:,i);
         g = g + betaa*G(:,i);
         w = w + betaa*W(:,i);
         for s = i+1:n-1
           betaa = -Q(:,s+1)'*z / c(s);
           z = z + betaa*D(:,s);
           g = g + betaa*G(:,s);
           w = w + betaa*W(:,s);
         end   
      end
 
      w = r + rhoo*w;
      betaa = -Q(:,1)'*w/ d;
      w = w + (rhoo*betaa)*W(:,n);
      g = g + w + betaa*G(:,n);

      for s = 1:i-1
        betaa = -Q(:,s+1)'*w / c(s);
        g = g + betaa*G(:,s);
        w = w + betaa*D(:,s);
      end

      D(:, i) = w - u;
      c(i) = Q(:,i+1)'*D(:,i);
      
      if c(i) == 0
          flag = -1;
          return
      end
      
      G(:,i) = g;
      g_tld = M\g;
      W(:,i) = A*g_tld;

      alphaa = f/c(i);
      if i < n - 1
         u = u - alphaa*D(:,i);
      end
      
      alphaa = rhoo*alphaa;
      x = x + alphaa*g_tld;
      r = r - alphaa*W(:,i);
 
      err = norm(r)/bnrm2;
      iter = iter + 1;
 
      if err < tol
          flag = 0;
          return
      end

      if iter >= max_it
         return
      end
    end
    
    e = Q(:,1)'*r;
    betaa = - e / d;
    w = r + (rhoo*betaa)*W(:,n);
    g = w + betaa*G(:, n);

    for s = 1 : n-1
       betaa = - Q(:,s+1)'*w / c(s);
       g = g + betaa*G(:,s);
       w = w + betaa*D(:, s);
    end

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

    if c(n) == 0
       flag = -1;
       return
    end
  end