function [x,it] = pcg_(A,b,prec_type,eeps,maxit,par)
n  = size(A,1);
x  = zeros(n,1);  % x0  then for x^{k+1}


if nargin < 6,   par   = 1;            end
if nargin < 5,   maxit = 100;          end
if nargin < 4,   eeps  = 1e-4;         end
if nargin < 3,   type  = 'none';       end

r0 = b-A*x;   
% z0 = Pinv * r0
switch lower(prec_type)
  case 'jacobi'
    z0 = (1./diag(A)).*r0;
  case 'cholinc'
    %L  = ichol(A,struct('type','ict','droptol',4e-1)); 
    L  = ichol(A);
    z0 = L'\(L\r0);
  case 'exact'
    z0 = A\r0;
  otherwise % 'none'
    z0 = r0;
end
r = r0;   z = z0;   rz_old = r'*z;
p = z0;
it = 0;

while (norm(r)/norm(r0) > eeps ) && (it < maxit)
  it = it + 1;
  Ap  = A*p;
  pAp = p'*Ap; 
  al = rz_old/pAp;
  x  = x + al*p;
  r  = r - al*Ap;
  switch prec_type % z = Pinv * r
    case 'jacobi'
      z = (1./diag(A)).*r;
    case 'cholinc'
      z = L'\(L\r);
    case 'exact'
      z = A\r;
    otherwise % 'none'
      z = r;
  end
  rz = r'*z;
  be = rz/rz_old; rz_old = rz;   % be = -(z'*Ap)/pAp;
  p  = z + be*p;
end

end