pcg.m

五点差分型多重网格方法：各种插值算子的比较）

%PCG    Preconditioned conjugate gradient method
%
%       [X,RESIDS,ITS]=PCG(A,B,X0,RTOL,PRTOL,MAX_IT,MAX_TIME,MAX_MFLOP)
%       solves the system AX = B using the preconditioned conjugate gradient 
%       method with the given tolerances and limits.  A should be symmetric
%       positive definite.
%

% James Bordner and Faisal Saied
% Department of Computer Science
% University of Illinois at Urbana-Champaign
% 10 April 1995

% Modified for Matlab Version 6 Compatability
% Ryan McKenzie
% University of Kentucky Center for Computational Sciences
% April 2004
%
% For some reason, locally generated variables cannot be seen outside the scope of a
% particular function in version 6 unless they have a global reference. This seems to 
% only occur when the newly generated variable is passed as a parameter. I have taken
% locally generated variables throughout MGLab and "bridged" them to their destination
% routines using global references. It's an ugly fix, so maybe someone should come up
% with a more centralized solution.
%
% Removed stopping criteria and other solver information from the parameter list since
% they are all stored in global constants anyway.

function [x,resids,its] = pcg(A,b,x0)

include_globals
include_bridge_globals

% gobally referencing variables from the solve routine "bridge"
b = b_in_sol_method;
A = A_in_sol_method;
x0 = x_in_sol_method;

x = x0; 
r = b - A*x; 
z = precondition(A,r);
prn = norm(z);

iter = 0;
results=update_results([],'PCG',iter,prn);

bn = 1;
rn = 1;
pbn = precondition(A,b);

% z could be named pr but we use the more common notation.

while (~converged(bn,pbn,rn,prn,iter,rtol,prtol,max_it,max_time))

   rho = r' * z;
   if iter == 0
      p = z;
   else
      beta = rho / rho_old;
      p = z + beta * p;
   end
   rho_old = rho;

   Ap = A*p;
   mu = p' * Ap;
   alpha = rho / mu;
   x = x + alpha *  p;
   r = r - alpha * Ap;
   z = precondition (A,r);

   prn = norm(z);

   iter = iter+1;
   results=update_results(results,'PCG',iter,prn);
end

its=results(:,1);
resids=results(:,4);