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