function [sol,iter,r0,Iterates] = mycgs (x,A,b,eps,kmax,cgflag)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MYCGS - Conjugate Gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 451-551
%% this function solves the system Ax=b using Conjugate Gradient method
%% refs: [C. T. Kelley, "Iterative Solution of Linear and Nonlinear Equations"]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2004/10
  %% changes added by Doug Fettig, 1/2005 for interface with cgdemo.m
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% INSTRUCTIONS: 
%% calling sequence: 
%%     [sol,iter,r0] = mycgs (x,A,b,eps,kmax)
%%% or
%%     [sol,iter,r0,I] = 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
%% 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.
%%%% optional parameters:
%% cgflag is a flag to switch to steepest descent method
%% if cgflag ==0 or no flag is set then the default (CG) method will be used
%% if cgflag != 0 then the steepest descent method is used, 
%%%%% 
%% sol is the solution 
%% iter is the actual number of iterations 
%% r0 on output is the norm of residual
%% Iterates contains all iterates 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% If A is not positive definite or not symmetric, then the routine
%%% outputs an error message and terminates.

spd = 1;
if isequal(A,A') 
  E = eig(A);
  for n = 1 : (size(E,1)) if ( E(n) <= 0 ) spd = 0;end; end;
else 
  spd = 0;
end;
if spd ==0             
  Error = [' INPUT ERROR: A IS NOT SYMMETRIC POSITIVE DEFINITE    '
           'Please use a symmetric, positive definite matrix for A']
  Iterates = 0; sol=0; iter=0; r0=0;
  return
end
    
%%% checks cgflag: if none, default is CG
if nargin < 6                          
    cgflag = 0;
end

%%% begin algorithm
                                        %%% initialize residual,its norm,k
r  = b - A*x;
r0 = norm(r)^2;
k  = 1; 

if norm(b) == 0         
    normb = eps;        
else 
    normb = eps * norm(b);      
end                                   

%%% initialize Iterates
Iterates = zeros(size(A,1), kmax);  Iterates (:,1) = x;

                                      %%% start the loop
while ( sqrt(r0) > normb && k < kmax) 

  if (k==1) p = r;
  else 
         beta = r0/rprev;               
                                      %% check if cg method or if steepest                           
         if cgflag ~= 0                     
             p = r;                   %% steepest descent
         else 
             p = r + beta*p;          %% new search direction for cg
         end  
     end;

  w = A * p;                          %% auxiliary vector 
  alpha = r0 / (p' * w);              %% found optimal alpha in line search
  x = x + alpha*p;                    %% new guess
  r = r - alpha*w;                    %% the residual is in fact r=b-A*x
  rprev = r0;                         %% ... save norm of the old residual    
  r0 = ( norm(r) )^2;                 %% ... compute norm of new residual
  k = k + 1;                          %% ... increase number of iterations
  Iterates(:, k) = x;                 %% store current iterate
end;

sol = x; iter = k; r0=sqrt(r0);       
                                     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
