nikcleju@0: % cgsolve.m
nikcleju@0: %
nikcleju@0: % Solve a symmetric positive definite system Ax = b via conjugate gradients.
nikcleju@0: %
nikcleju@0: % Usage: [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)
nikcleju@0: %
nikcleju@0: % A - Either an NxN matrix, or a function handle.
nikcleju@0: %
nikcleju@0: % b - N vector
nikcleju@0: %
nikcleju@0: % tol - Desired precision.  Algorithm terminates when 
nikcleju@0: %    norm(Ax-b)/norm(b) < tol .
nikcleju@0: %
nikcleju@0: % maxiter - Maximum number of iterations.
nikcleju@0: %
nikcleju@0: % verbose - If 0, do not print out progress messages.
nikcleju@0: %    If and integer greater than 0, print out progress every 'verbose' iters.
nikcleju@0: %
nikcleju@0: % Written by: Justin Romberg, Caltech
nikcleju@0: % Email: jrom@acm.caltech.edu
nikcleju@0: % Created: October 2005
nikcleju@0: %
nikcleju@0: 
nikcleju@0: function [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)
nikcleju@0: 
nikcleju@0: if (nargin < 5), verbose = 1; end
nikcleju@0: 
nikcleju@0: implicit = isa(A,'function_handle');
nikcleju@0: 
nikcleju@0: x = zeros(length(b),1);
nikcleju@0: r = b;
nikcleju@0: d = r;
nikcleju@0: delta = r'*r;
nikcleju@0: delta0 = b'*b;
nikcleju@0: numiter = 0;
nikcleju@0: bestx = x;
nikcleju@0: bestres = sqrt(delta/delta0); 
nikcleju@0: while ((numiter < maxiter) && (delta > tol^2*delta0))
nikcleju@0: 
nikcleju@0:   % q = A*d
nikcleju@0:   if (implicit), q = A(d);  else  q = A*d;  end
nikcleju@0:  
nikcleju@0:   alpha = delta/(d'*q);
nikcleju@0:   x = x + alpha*d;
nikcleju@0:   
nikcleju@0:   if (mod(numiter+1,50) == 0)
nikcleju@0:     % r = b - Aux*x
nikcleju@0:     if (implicit), r = b - A(x);  else  r = b - A*x;  end
nikcleju@0:   else
nikcleju@0:     r = r - alpha*q;
nikcleju@0:   end
nikcleju@0:   
nikcleju@0:   deltaold = delta;
nikcleju@0:   delta = r'*r;
nikcleju@0:   beta = delta/deltaold;
nikcleju@0:   d = r + beta*d;
nikcleju@0:   numiter = numiter + 1;
nikcleju@0:   if (sqrt(delta/delta0) < bestres)
nikcleju@0:     bestx = x;
nikcleju@0:     bestres = sqrt(delta/delta0);
nikcleju@0:   end    
nikcleju@0:   
nikcleju@0:   if ((verbose) && (mod(numiter,verbose)==0))
nikcleju@0:     disp(sprintf('cg: Iter = %d, Best residual = %8.3e, Current residual = %8.3e', ...
nikcleju@0:       numiter, bestres, sqrt(delta/delta0)));
nikcleju@0:   end
nikcleju@0:   
nikcleju@0: end
nikcleju@0: 
nikcleju@0: if (verbose)
nikcleju@0:   disp(sprintf('cg: Iterations = %d, best residual = %14.8e', numiter, bestres));
nikcleju@0: end
nikcleju@0: x = bestx;
nikcleju@0: res = bestres;
nikcleju@0: iter = numiter;
nikcleju@0: