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