%----------------------------------------------
%     The conjugate gradient method (cgm).
%            [A]{x} - {f}  =  0
%----------------------------------------------
function  x = CG(A,x,f,neq)

%     --------------------------
%     Calculate initial residual
%         {r} = {f} - [A]{x}
%     --------------------------
      r = f -A*x;

%     ----------------------------------------
%     Set:   
%        ep   = convergence test value
%        iter = iteration counter
%        rmag = initial magnitude of residual
%     ----------------------------------------
      ep   = 1.0e-05;         
      rmag = r'*r;
      iter = 0;

%     ---------------------
%     Begin iterations
%     ---------------------
      while  rmag > ep & iter <= neq+1 
         iter = iter +1;
%        ------------------------------------
%        Determine new search direction {p}
%        ------------------------------------
         if iter == 1;
%           ---------------------------------
%	    Use steepest decent, {r}
%	    {p} = {r}
%           ---------------------------------
            p = r;
	 else
%           ----------------------------------------
%	    Make new direction conjugate to last {p}
%	    Beta = |r|{Ap}/pAp
%	    {p} = {r} - Beta*{p}
%           ----------------------------------------
            rAp  = r'*Ap;
            Beta = rAp/pAp;
            p    = r - Beta*p;
         end

%        -----------------------------------------------------
%        Determine distance along {p} to minimum point
%        rp = |r|{p};  {Ap} = [A]{p};  pAp = |p|{Ap}
%        -----------------------------------------------------
         rp = r'*p;
	 Ap = A*p;
	 pAp = p'*Ap;
         s = (rp)/(pAp);
	 x = x + s*p;
        
%        ------------------------------------------------------
%	 Determine new residual (direction of steepest decent)
%                {r} = {f} - [A]{x}
%        or, update old residual direction (faster)
%                {r} = {r} - s*{Ap}
%        ------------------------------------------------------
	 r = r -s*Ap;
%        --------------------------------
%        Determine magnitude of residual
%        rmag = |r|{r}
%        --------------------------------
         rmag = r'*r;
      end

      if iter <= neq
            fprintf(1,'\n ')
            fprintf(1,'\n     --------------------------------------')
            fprintf(1,'\n     Converged     ')
            fprintf(1,'\n     rmag,iter = %10.3e,%5d',rmag,iter) 
            fprintf(1,'\n     --------------------------------------')
            fprintf(1,'\n ')
	 else
            fprintf(1,'\n ')
            fprintf(1,'\n     --------------------------------------')
            fprintf(1,'\n     Did not converge     ')
            fprintf(1,'\n     rmag,iter = %10.3e,%5d',rmag,iter) 
            fprintf(1,'\n     --------------------------------------')
            fprintf(1,'\n ')
	 end

