function [x, error, total_iters, it_hist] = ...
                     pcgsol(x0, b, atv, params, pcv)
% Preconditioned Conjugate-Gradient solver
%
% C. T. Kelley, December 12, 1993
%
% This code comes with no guarantee or warranty of any kind.
%
% function [x, error, total_iters, it_hist] 
%                    = pcgsol(x0, b, atv, params, pcv)
%
%
% Input:        x0=initial iterate
%               b=right hand side
%               atv, a matrix-vector product routine
%			atv must return Ax when x is input
%			the format for atv is
%                       function ax = atv(x)
%               params = two dimensional vector to control iteration
%                        params(1) = relative residual reduction factor
%                        params(2) = max number of iterations
%               pcv, a routine to apply the preconditioner
%                       if omitted, the identity is used.
%                       The format for pcv is
%                       function px = pcv(x).
%
% Output:       x=solution
%               error = vector of iteration residual norms
%               total_iters = number of iterations
%               it_hist (optional) = matrix of all iterations
%			useful for movies
%
% 
%

%
% initialization
%
if nargout == 4 it_hist=[]; end
n=length(b); errtol = params(1); maxiters = params(2); error=[]; x=x0;
if nargout == 4; it_hist=[it_hist, x]; end
r=b - feval(atv, x);
if nargin == 4
    z=r;
else
    z = feval(pcv, r);
end
rho=z'*r;
tst=norm(r);
terminate=errtol*norm(b);
error=[error,tst];
it=1;
while((tst > terminate) & (it <= maxiters))
%
%
%
if(it==1) 
	p = z;
else
	beta=rho/rhoold;
	p = z + beta*p;
%
% end if
%
end
w = feval(atv, p);
alpha=p'*w;
%
% Test here to make sure the linear transformation is positive definite.
% A non-positive value of alpha is a very bad sign.
%
if(alpha <= 0)
    [alpha, rho, it]
    error(' negative curvature ') 
end
alpha=rho/alpha;
x=x+alpha*p;
if nargout == 4; it_hist=[it_hist, x]; end
r = r - alpha*w;
tst=norm(r);
rhoold=rho;
if nargin == 4
    z=r;
else
    z = feval(pcv, r);
end
rho=z'*r;
it=it+1;
error=[error,tst];
%
% end while
%
total_iters=it-1;
end