
function [x, res, iter, flag] = mycg( A, b, M, varargin ) 
%
%  MYCG:
%     computes an approximate solution of A x = b using the 
%     conjugate gradient algorithm.
%
%  Usage:
%
%  [x, res, iter, flag] = mycg( A, b, M, ) 
%  [x, res, iter, flag] = mycg( A, b, M, x ) 
%  [x, res, iter, flag] = mycg( A, b, M, x, options )   
%
%
%  Input parameters
%
%  A       the system matrix A.
%          is either a real(n,n) symmetric positive definite matrix 
%          the name of a function that computes y = A * x.
%          The calling sequence of this function is [y] = Ax(x).
%
%  b       real(n)
%          the right-hand side. 
%
%  M       the preconditioner
%          is either a symmetric positive definite n by n matrix or the 
%          name of a function that solves M z = x.
%          The calling sequence of this function is  [z] = M(x).
%          No preconditioning is applied when M = [].
%
%  x       real(n)
%          the initial approximate solution. Default x = 0.
%
%  options algorithm options
%          options.tol    stopping tolerance (Default tol = 1.e-7)
%                         the algorithm attempts to generate x such that   
%                         || M^(-1) ( b-Ax ) || <= tol
%          options.maxit   maximum number of iterations (Default maxit =
%          100)
%
%
%  Output parameters: 
%
%  x       Final approximate solution. 
%
%  flag    integer
%          Error flag, as follows: 
%              0 => normal termination; desired tolerance met. 
%              1 => desired tolerance not met after maxit iterations.
%              2 => in iteration iter an inner product of the form
%                   v' M^(-1) v was not positive, so the preconditioning 
%                   matrix M  does not appear to be positive definite.
%                   x  will contain the last CG iterate computed before
%                   this event. 
%              3 => in iteration iter an inner product of the form
%                   q' A q was not positive, so the matrix A  does not 
%                   appear to be positive definite.
%                   x   will contain the last CG iterate computed before
%                   this event. 
%
%  res     real
%          vector of residual norms. 
%          res(i+1) is the residual norm of iteration i.
%
%  iter    integer
%          contains the total number of iterations in CG.
%          the number of  matrix vector products is given by  iter + 1.
% 
%
%
%
%  Matthias Heinkenschloss  
%  Feb 13, 2009
%   
%

%% set algorithm parameters
if( nargin < 4 );
    x = zeros(size(b,1),1);
else
    x = varargin{1};
end
if( nargin < 5 );
    options = [];
else
    options = varargin{2};
end
maxit   = 100;
tol     = 1.e-7;
if( ~isempty(options) )
    if( isfield(options,'maxit') );    maxit   = options.maxit; end
    if( isfield(options,'tol') );      tol     = options.tol; end
end

%% initialization
flag = 0;
iter = 0;

% compute initial residual
if( isstr(A) )
    r = b - feval( A, x );
else
    r = b - A*x;
end

%% begin iteration
for iter = 0:maxit
    
    % apply preconditioner
    if( isstr(M) )
        z = feval( M, r );
    elseif( ~isempty(M) )
        z = M\r;
    else
        z = r;
    end
    rho = (r'*z);
    if( rho < 0 ); flag = 2; break, end
    res(iter+1) = sqrt(rho);
    % check convergence
    if ( res(iter+1) <= tol ), return, end
    
    % direction vector
    if ( iter > 0 ),
        beta = rho / rho_1;
        p = z + beta*p;
    else
        p = z;
    end
    
    if( isstr(A) )
        q = feval( A, p );
    else
        q = A*p;
    end
    alpha = (p(:)'*q(:));
    
    if( alpha <= 0 ); flag = 3; break, end
    alpha = rho / alpha;
    
    % update approximation vector
    x = x + alpha * p;
    
    % compute residual
    r = r - alpha * q;
    
    rho_1 = rho;
    
end


% no convergence
flag = 1;

% END mycg.m
