cg_for_IP.m

function [x, error, iter, flag] = cg_for_IP(funAx, x, b, funPrec, max_it, tol)

% Updated by MK for POEMA and Interior Point algorithm:
%
% funAx is a function for matrix-vector multiplication
% funPrec is a function that solves the system Mz = r in preconditioning
%

%  -- Iterative template routine --
%     Univ. of Tennessee and Oak Ridge National Laboratory
%     October 1, 1993
%     Details of this algorithm are described in "Templates for the
%     Solution of Linear Systems: Building Blocks for Iterative
%     Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
%     Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
%     1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
%
%  [x, error, iter, flag] = cg(A, x, b, M, max_it, tol)
%
% cg.m solves the symmetric positive definite linear system Ax=b 
% using the Conjugate Gradient method with preconditioning.
%
% input   A        REAL symmetric positive definite matrix
%         x        REAL initial guess vector
%         b        REAL right hand side vector
%         M        REAL preconditioner matrix
%         max_it   INTEGER maximum number of iterations
%         tol      REAL error tolerance
%
% output  x        REAL solution vector
%         error    REAL error norm
%         iter     INTEGER number of iterations performed
%         flag     INTEGER: 0 = solution found to tolerance
%                           1 = no convergence given max_it

  flag = 0;                                 % initialization
  iter = 0;

  bnrm2 = norm( b );
  if  ( bnrm2 == 0.0 )
      bnrm2 = 1.0;
  end
  
  x = sparse(x);

  r = b - funAx(x);
  error = norm( r ) / bnrm2;
  if ( error < tol )
      return
  end

  for iter = 1:max_it                       % begin iteration

     %z  = M \ r;
     z = funPrec(r);
     rho = (r'*z);

     if ( iter > 1 )                       % direction vector
        beta = rho / rho_1;
        p = z + beta*p;
     else
        p = z;
     end

     q = funAx(p);
     alpha = rho / (p'*q );
     x = x + alpha * p;                    % update approximation vector

     r = r - alpha*q;                      % compute residual
     error = norm( r ) / bnrm2;            % check convergence
     if ( error <= tol )
         break
     end 

     rho_1 = rho;

  end

  if ( error > tol )
      flag = 1;
  end         % no convergence

% END cg.m