function [sol,iter,rnorm] = mycg (x,A,b,tol,kmax)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MYCGS - Conjugate Gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 654
%% refs: [C. T. Kelley, "Iterative Solution of Linear and Nonlinear Equations"]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2009/10
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% INSTRUCTIONS: 
%% calling sequence: 
%%     [sol,iter,r0] = mycgs (x,A,b,eps,kmax)
%% x is the initial guess (a vector)
%% A is the system matrix (positive definite, symmetric)
%% b is the right hand side
%% myeps, 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 
%% iter is the actual number of iterations 
%% rnorm on output is the norm of residual
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r  = b - A*x;                       %% SGEMV followed by SAXPY
r0 = r'*r;                          %% SDOT or SNRM2
rnorm = sqrt(r0);

k  = 1; 
p = r;                                %% SCOPY or SAXPY
                               
while ( rnorm > tol && k < kmax) 
  w = A * p;                          %% SGEMV or SSBMV 
  aconj = p' * w ;                    %% SDOT
  
  
  if norm(aconj) < 1e-8, break; end   %% SNRM2 ?
  alpha = r0 / aconj;                 %% optimal alpha in line search
  x = x + alpha*p;                    %% SAXPY
  r = r - alpha*w;                    %% SAXPY
  
  rprev = r0;                            
  r0 = r'*r;                          %% SNRM2 or SDOT
  rnorm = sqrt(r0);   
 
  beta = r0/rprev;                                         
  p = r + beta*p;                     %% update search direction. SAXPY
  k = k + 1;
  plot(x); hold on;
end;

sol = x; iter = k;    
                  
pause,hold off;clf;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
