function [sol,iter,r0] = mypcg (x,A,b,M,eps,kmax)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MYCGS - Conjugate Gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 451-551
%% this function solves the system Ax=b using Preconditioned 
%% Conjugate Gradient method, see mycgs.m
%% refs: [C. T. Kelley, "Iterative Solution of Linear and Nonlinear Equations"]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2004/10
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% INSTRUCTIONS: 
%% calling sequence: 
%%     [sol,iter,r0] = mypcgs (x,A,b,M,eps,kmax)
%% x is the initial guess (a vector)
%% A is a positive definite, symmetric matrix
%% b is the right hand side
%% M is the preconditioner (nonsingular matrix "close" to A
%% eps, kmax determine the stopping criterium: 
%%    algorithm stops when the norm
%%    of current residual is less than eps times the norm of b on input
%%    or if the number of iterations exceeds kmax 
%% on output, sol is the solution 
%% k is the actual number of iterations 
%% r0 on output is the norm of residual
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                                      %%% initialize residual,its norm,k
r = b - A*x;
z = M \ r;

r0 = norm(r)^2;
k  = 1; 
normb = eps*norm(r);                  %%% precompute the bound of residual size

                                      %%% start the loop
while ( sqrt(r0) > normb && k < kmax) 
  z = M \ r;                          %% action of preconditioner 
  tau = z'*r;                         %% ... compute new coefficient

  if (k==1) p = z;
  else 
         beta = tau/tauprev; 
         p = z + beta*p;              %% new search direction
  end;

  w = A*p;                            %% auxiliary vector 
  alpha = tau/(p'*w);                 %% found optimal alpha in line search
  x = x + alpha*p;                    %% new guess
  r = r - alpha*w;                    %% recompute residual


  tauprev = tau;                      %% ... save previous coefficient
  r0   = r'*r;                        %% ... compute norm of new residual
  k = k+1;                            %% ... increase number of iterations

end;

sol = x; iter = k; r0=sqrt(r0);       %% on output, sol is the solution 
                                      %% k is the actual number of iterations 
                                      %% r0 on output is the norm of residual

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

