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