function [x, error, iter, flag] = pbicgstab(A, x, b, M, max_it, tol, M1, M2)
% solves the linear system Ax=b using the
% BiConjugate Gradient Stabilized Method with preconditioning.

  iter = 0;% initialization
  flag = 0;

  bnrm2 = norm( b );
  if  ( bnrm2 == 0.0 ), bnrm2 = 1.0; end

  r = b - A*x;
  error = norm( r ) / bnrm2;
  errorhist = [ ];                   
  errorhist(1) = error;              

  if ( error < tol ) return, end

  omega  = 1.0;
  r_tld = r;

  for iter = 1:max_it,                              

     rho   = ( r_tld'*r );                          
     if ( rho == 0.0 ) break, end

     if ( iter > 1 ),
        beta  = ( rho/rho_1 )*( alp/omega );
        p = r + beta*( p - omega*v );
     else
        p = r;
     end
     p_hat = M2 \ (M1 \ p);    
     v = A*p_hat;
     alp = rho / ( r_tld'*v );
     s = r - alp*v;
     if ( norm(s) < tol ),                          
        x = x + alp*p_hat;
        resid = norm( s ) / bnrm2;
        break;
     end
     
     s_hat = M2 \ (M1 \ s) ;     
     t = A*s_hat;
     omega = ( t'*s) / ( t'*t );

     x = x + alp*p_hat + omega*s_hat;            

     r = s - omega*t;
     error = norm( r ) / bnrm2;                    
     errorhist(iter+1) = error;            

     if ( error <= tol ), break, end

     if ( omega == 0.0 ), break, end
     rho_1 = rho;

  end

  if ( error <= tol | s <= tol ),                   
     if ( s <= tol ),
        error = norm(s) / bnrm2;
        errorhist(iter+1) = error;         
     end
     flag =  0;
  elseif ( omega == 0.0 ),                         
     flag = -2;
  elseif ( rho == 0.0 ),
     flag = -1;
  else                                             
     flag = 1;
  end

  error = errorhist;  